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Preface 



Manual Scope 



Audience 



What is in This Manual 



Sun Microsystems provides a variety of software and hardware floating-point 
options. Some options provide fast performance; most conform substantially, 
and all conform at least partially, to the requirements of the ANSI/IEEE Std 
754-1985, the IEEE Standard for Binary Floating-Point Arithmetic. 

This manual describes features of Sun Microsystem’s various available floating- 
point implementations, including details of their use, performance, and confor- 
mance to the IEEE floating-point arithmetic standard. 

Note: 

All specific claims for conformance, accuracy, and speed are based on tests 
and measurements made using a preliminary version of Sun Software 
Release 3.2. Results obtained with the final Release 3.2 or with other 
releases may vary. 

This manual is limited to information not easily obtainable elsewhere; you 
should be familiar with the ANSI/IEEE Standard 754-1985, with the MC6888 1 
floating-point processor (if you have a Sun-3), and the Weitek WTL 1 164/1 165 
chip set (if you have a Floating-Point Accelerator or FPA). 

If you want or need more information than is contained in this manual, see one or 
more of the sources listed in the bibliography later in this preface. 

The description of floating-point options on Sun Microsystems computers is 
divided into four chapters and eight appendices: 

Chapter 1 describes the hardware and software options that can affect floating- 
point numeric operations on Sun Microsystems computers including hardware, 
compiler code generation, and expression evaluation. 

Chapter 2 describes floating-point numerics as applied to Sun Micosystems 
hardware and software. This includes discussions of rounding errors and differ- 
ring numerical results depending on different floating-point options, using 
floating-point operations from compiled and assembly language programs, and 
operational details of Sun Microsystems’ floating-point implementations. 

Chapter 3 describes aspects of the ANSI/IEEE Standard 754-1985, and Sun 
Microsystems floating-point conformance to the standard. 
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References for Further 
Information 



Chapter 4 describes benchmarks that have been and can be used to test the per- 
formance and conformance of floating-point numerics on Sun Microsystems 
computers. 

Appendix A describes changes made to adb to support the FPA, examples of 
FPA disassembly, and FPA register use. 

Appendix B describes changes made to dbx and dbxtool to support the 
FPA, examples of FPA disassembly, and FPA register use. 

Appendix C describes the FPA assembly-language syntax supported by as in 
Release 3.1 and later. 

Appendix D describes a number of IEEE floating-point functions that are not 
usually available in higher-level languages. 

Appendix E describes the input files used by the SPICE benchmark program. 

Appendix F describes the differences in function between the two different mask 
versions of the MC6888 1 coprocessor used by Sun’s floating-point processors. 

Appendix G describes assembly-level inline code expansions that let you 
integrate assembly-language routines into C, FORTRAN, and Pascal programs. 

Appendix H describes aspects of the floating-point implementation affected by 
compliance with the UNIX System V interface. 

The following provide more information about Sun’s language processors: 

FORTRAN Programmer s Guide, 800-1371-02, Sun Microsystems, 1986. 

Pascal Programmer’ s Guide, 800-1376-01, Sun Microsystems, 1986. 

Assembly Language Reference Manual, 800-1372-02, Sun Microsystems, 
1986. 

The following provide more information about Sun’s floating-point hardware: 

Sun Floating-Point Accelerator User’ s Manual, 800-1378-02, 1986. 

Preliminary Data, WTL 1164/1165 Low-Latency 64-bit IEEE Floating Point 
Multiplier! ALU, 231669-001, Weitek Corp., Sunnyvale, CA, 1986. 

MC68881 Floating-Point Coprocessor User’ s Manual, MC6888 1UM/ AD, 
Prentice-Hall, 1985. 

MC68020 32-Bit Microprocessor User’ s Manual, second edition, 
MC68020UM/AD REV 1, Motorola Inc., 1985. 

Fast Floating-Point Processor Integration Manual, Sky Computers, Lowell, 
MA, October 1984. 

The following provide more information about the IEEE Standard: 

IEEE Standard for Binary Floating-Point Arithmetic, ANSI/IEEE Std 754- 
1985, IEEE, New York, 1985. 

Apple Numerics Manual, Addison-Wesley, 1986. 



- xii - 




Preface — Continued 



Test Program Information 



Coonen, Contributions to a Proposed Standard for Binary Floating-Point 
Arithmetic, Ph. D. thesis. University of California, Berkeley, 1984. 

Cody et al., "A Proposed Radix- and Word-length-independent Standard for 
Floating-Point Arithmetic," IEEE Computer, August 1984. 

Stevenson et al., Cody, Hough, Coonen, various papers proposing and 
analyzing a draft standard for binary floating-point arithmetic, IEEE Com- 
puter, March 1981. 

Coonen, "An Implementation Guide to a Proposed Standard for Floating- 
Point Arithmetic," IEEE Computer, January 1980. 

The Proposed IEEE Floating-Point Standard, special issue of the ACM SIG- 
NUM Newsletter, October 1979. 

The following defines the UNIXt System V Interface: 

System V Interface Definition, Issue 2, AT&T, Indianapolis, IN, 1986. 

The following provide more information about testing and error analysis: 

Dongarra, Performance of Various Computers Using Standard Linear Equa- 
tions Software in a FORTRAN Environment, Argonne National Laboratory, 
Argonne, IL, 1986. 

Karpinski, "Paranoia: a Floating-Point Benchmark," Byte, February 1985. 

Spafford, Flaspohler, A Report on the Accuracy of some Floating Point Math 
Functions on Selected Computers, Technical Report GIT-CS 85/06, Georgia 
Institute of Technology, Atlanta, 1985. 

"Appendix: Accuracy of Numerical Calculations," in HP-15C Advanced 
Functions Handbook, 00015-90011, Hewlett-Packard, 1982. 

Cody and Waite, Software Manual for the Elementary Functions, Prentice- 
Hall, 1980. 

Bunch, Dongarra, Moler, Stewart, Linpack Users’ Guide, SIAM, Philadel- 
phia, 1979. 

Cumow, Wichmann, "A Synthetic Benchmark," The Computer Journal, 
British Computer Society, London, February 1976. 

Sterbenz, Floating-Point Computation, Prentice-Hall, 1974. 

For information on machine-readable source of the IEEE Test Vectors or the 
SPICE program, contact Cindy Manley, EECS/ERL Industrial Support Office, 
461 Cory Hall, University of California, Berkeley, CA 94720. 

For information on machine-readable source for the Paranoia test program, 
contact Richard Karpinski at ucbvaxlucsfcgl! cca.ucsf! dick. 

For information on Alex Liu’s test programs, contact him at 
zliu@weyl.berkeley.edu. 



t UNIX is a trademark of AT&T Bell Laboratories. 
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Netlib provides a numerical source code distribution service from which many 
popular benchmark programs can be obtained, including paranoia, 
elefunt, Unpack, and whetstone, as well as the current version of 
Dongarra’s Argonne Technical Memorandum containing Linpack benchmark 
results. For netlib instructions, send a message "send netlib-paper from misc" to: 
research! netlib or netlib@anl-mcs.arpa. 

For information on the FCVS tests, write to Federal Software Testing Center, 
Suite 1 100, 5203 Leesburg Pike, Falls Church, VA 22041. 
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Sun Floating-Point Options 



1.1. Hardware Floating- 
Point Options 

Sun-2 Systems 



Sun-3 Systems 



1.2. Compiler Code 
Generation Options 



This manual describes floating-point support for Sun Microsystems’ Sun-2 and 
Sun-3 systems. Sun’s earlier models are not discussed here: The phrase "runs on 
any Sun" means that programs compiled fora Sun-2 ran on any Sun-2 or Sun-3, 
while programs compiled for a Sun-3 run on any Sun-3. Throughout this 
manual, mention of "all" Sun systems refers to Sun-2’s and Sun-3’s. 

The floating-point hardware options usable with a Sun Microsystems computer 
depend on the model Sun being used, either a Sun-2 or Sun-3. 

Sun-2 systems use a 10 MHz MC68010 as their CPU and use either a Multibus 
or VME system bus. The Sky Fast Floating-Point Processor (FFP) board exists 
in both Multibus and VME versions for the Sun-2. Programs compiled for Sun- 
2’s will usually ran on Sun-3’s as long as they do not exploit any specific 
features unique to Sun-2’ s, such as the Sky FFP. 

Sun-3 systems use either a 15 or 16.7 MHz MC68020 as their CPU and use a 
VME system bus. An MC6888 1 floating-point coprocessor is standard on some 
models and optional on others. Sun-3’s use MC6888 l’s of either the mask set 
A79J, running at 12.5 MHz, or of the mask set A93N, running at 16.7 MHz. 

Also available on some Sun-3 models is the Sun Floating-Point Accelerator 
(FPA), which is an optional board based upon the Weitek 1 164/1165 chip set. 

The FPA requires the presence of an MC6888 1. 

Programs compiled for Sun-3’s will not ran on Sun-2’s. 

Sun’s compilers for FORTRAN (f 7 7), C (cc), and Pascal (pc) use the options 
described below in the same way. With each of these compilers, if you compile 
and link programs in separate steps, you should usually specify the same set of 

options at each step. For example, for FORTRAN: 
— 

f77 -pg -0 -ffpa -c parti. f 
f77 -pg -O -ffpa -c part2.f 

f77 -pg -0 -ffpa partl.o part2.o -o executable . out 
v. > 



Note: 

None of Sun’s compilers attempt extensive code optimization unless 
requested. 
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Full Optimization (-0) 



Partial Optimization (-P) 



Optimization Side Effects 



Inline Expansion 



The -0 option is the usual way to request maximum code optimization. In 
Release 3.2 f77, -O invokes /usr/lib/iropt to perform the following 
global optimizations on the program: 

□ Invariant code removal from loops 

□ Induction variable strength reduction 

□ Common subexpression elimination 

□ Copy propagation 

□ Register allocation 

□ Dead code elimination 

Note that iropt’s global optimizations sometimes cause conformance prob- 
lems with the IEEE floating-point standard. 

Floating-point variable register allocation is discussed below under the 
-f store option. 

In Release 3.2 cc and pc, -0 invokes a more limited set of global optimiza- 
tions. In all three compilers, -0 also invokes / lib/ c2 to perform local 
optimizations on the compiled code before assembly. 

Partial optimization (-P) generates smaller data structures but often produces 
code almost as good as -0. 

Partial optimization is called for in cases where -O results in inordinately long 
compilation times: some modules generate huge internal data structures in the 
global optimizer ir opt that can lead to enormous compilation times with -0. 
If you think that is happening you can investigate this by printing out the 
compiler’s phases with the -v option; if compilation seems to hang up in 
iropt, use -P instead of -0. 

Don’t use -P as a general substitute for -O since -O usually produces more 
efficient code. 

iropt’s global optimizations sometimes cause problems with conformance to 
the IEEE floating-point standard. Floating-point variable allocation to registers 
is discussed with the -fstore option. Dead code elimination, common subex- 
pression elimination, and copy propagation may cause source code statements to 
not be executed when expected, and consequently IEEE exceptions may not be 
raised or may be raised a different number of times or in different places than the 
source code suggests. 

For examples, see the subroutine next d in the section "How to Use Special 
Features of the IEEE Standard" in Chapter 3, and the results of the Paranoia 
benchmark for the MC6888 1 in Chapter 4. 

Inline expansion of external calls is an optional third phase of optimizatioa You 
can select it independently of other optimization choices by specifying one or 
more . il files along with the source files in a compilation. 

The inline subroutine expansion program / usr/lib/inline replaces sub- 
procedure calls with equivalent code inline in the calling procedure, using 



#sun 

V microsystems 



Revision A of 12 September 1986 




Chapter 1 — Sun Floating-Point Options 5 



Floating-Point Code 
Generation Options 



Default Floating-Point Code 
Generation 



templates contained in the . il files. User-defined . il files may be used in 
addition to or in place of Sun-provided standard . il files. See Appendix G for 
a discussion on using inline subroutine expansion. 

Sun supports two approaches to floating-point code generation: switched and in- 
line. Switched code determines at runtime what the fastest available floating- 
point hardware is and selects from among available library routines to exploit 
that hardware. Thus switched code runs on any Sun, with any or no floating- 
point hardware, and takes advantage of whatever hardware is present. The 
results of floating-point computations may vary depending on the hardware that 
was used. Select switched floating-point code generation at compile time with 
the compiler option -f switch. 

Select software floating-point code generation with the compiler option 
-f sof t . Code so compiled runs on any Sun and always produces the same 
numerical result. Code compiled with -f sof t contains calls to floating-point 
libraries and is about as fast as switched floating point on machines without float- 
ing point hardware. 

In-line code generation puts floating-point instructions in line with CPU instruc- 
tions; such code generated in-line for specific hardware, by avoiding the over- 
head of library calls, always runs faster than switched code running on the same 
hardware. Remember that in-line code only runs on machines with the specific 
hardware installed for which the code was compiled. 

The available floating-point code generation hardware options are: 

□ -fsky Requires Sky FFP on Sun-2. 

□ -f 68881 Requires MC6888 1 on Sun-3. 

□ -ffpa Requires MC6888 1 and Sun FPA on Sun-3. 

-f soft is the default floating-Point code generation mode, in the absence of 
other specifications. This was chosen so that programs that use floating point in 
a casual way will obtain the same numerical results on all Suns. 

You may change your default by defining the environment variable 

F LOAT_OP T I ON with a . cshrc or . login command like 


setenv FLOAT_OPTION f 68 8 81 

V J 



Mixing Floating-Point Code 
Generation Styles 



Other legal choices are f switch, fsoft, fsky, or ffpa. 

Bear in mind when converting makefiles and shell scripts, that Sun releases 
before 3.0 made switched floating-Point the default, with -fsky the only 
option. 

An executable program composed of several independently compiled modules 
can contain only one kind of in-line hardware floating-point. If this restriction is 
violated, one of the following error messages 
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r 


Undefined: f sky_used 


\ 




Undefined : f 6 8 8 8 l_used 






Undefined: ffpa_used 




V 




j 



will be issued by the compiler when it tries to link the modules together. 

Although the practice is not recommended, you can link together modules com- 
piled with -f sof t, -f switch, and any one type of in-line hardware floating 
point. 

The compilers defeat attempts to specify -f sky on a Sun-3 or -f f pa on a 
Sun-2. 



Libraries and Executable File 
Sizes 



For small and moderate sized programs, inline code generation creates the smal- 
lest executable files because the fewest library subroutines are required. - 
f sof t requires library subroutines for all operations, while switched floating- 
point (-f switch) requires library subroutines for software floating-point and 
all the inline options and so is the largest of all. There are, however, some code 
generation technicalities that affect these generalizations: 



Each compiler searches certain libraries by default: 



r 




A 






Libraries Searched 


Compiler 


Normal 


Profiling (-p or -pg) 


cc 


c 


CJ? 


til 


F77 111 


U77 m c F77_p I77_p U77j> m_jo c_p 


pc 


pc m c 


pc_p m_j? c_p 






j 



In this table m, for instance, is shorthand for option -lm or the equivalent file 
name /usr/lib/libm.a,exceptfor c, which corresponds to -lc or 
/usr/lib/libc . a. These libraries contain subroutines which are in turn 
called by code generated by the compilers. The compilers may generate code to 
call different subroutines, depending on the -f . . . option in effect at compiler 
time. 

In general, 1 ibm . a and 1 ib c . a contain different subroutines corresponding 
toeach -f... option, libpc.a, libF77.a, libI77.a,and 
libU77 . a do not; the subroutines in these libraries are therefore created with 
-f switch, and so if they are called, then much of switched floating point is 
linked into the final executable file. 



This means that the executable file may be larger than expected. In some cases 
the switched floating point may be bypassed by using the inline subroutine 
expansion files described in Appendix G. 

FORTRAN source-language constructs that generate calls to subroutines in 
libF77 .a include, for real or double precision, x**y, mod, atan2, cabs, 
and nint ; for complex or doublecomplex, almost all operators and functions. 
With the Sun-provided inline subroutine expansion files, switched floating-point 
subroutines are not linked in except for for complex and doublecomplex x/y, 
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Summary: Which Option 
When? 



1.3. Expression Evaluation 
Options 



FORTRAN Style 



Greatest Available Precision 
Style 



-f store 



x**y, sqrt, exp, log, sin, and cos. 

Compile programs with -f sof t that use floating point in a casual way or not at 
all; this is the usual default anyway. Compile programs on a Sun-2 with 
-f switch that use floating point intensively but that are intended to be run on 
all possible Sun-2 and Sun-3 hardware configurations. Compile programs that 
are intended to use particular floating-point hardware intensively with in-line 
code for that hardware. 

Remember that the Sun FPA can only support 32 simultaneous processes. Con- 
sequently, on systems with a Sun FPA, -f f pa and -f swit ch should not be 
used for programs that use floating point casually, reserving floating-point con- 
texts for intensive floating-point programs. 

Programs that exploit certain details of the IEEE Standard may require 
-fsoft, -f 68881, or -ff pa; see Chapter 3. 

Some languages prescribe precise evaluation rules for expressions like these: 


a + b 

c * (d + e) 

J 



Others do not. Different kinds of floating-point hardware also lend themselves to 
different methods. 

The traditional method of expression evaluation in FORTRAN is that operations 
are performed in the precision of the most precise operand. Thus, a + b is 
evaluated in single precision if a and b are both single-precision, and in double 
precision if either a or b is double-precision. 

Since floating-point expressions usually suffer rounding errors when evaluated, 
these errors may be minimized if expressions are evaluated in the highest avail- 
able precision. Thus many C compilers evaluate a + b in double precision 
even if a and b are both single-precision. 

Sun’s FORTRAN and Pascal compilers generally perform FORTRAN style expres- 
sion evaluation. An exception is when the -f68881 option is specified; then 
expressions are evaluated in the extended-precision registers of the MC6888 1. 
Sun’s C compiler generally performs double-precision expression evaluation; if 
-f 68881 is specified, evaluation is in extended precision. The -f store, 

-f single, and -f single2 options, described below, provide alternate 
methods. 

When the -f 68881 option is used in conjunction with -0 or -P, the com- 
pilers allocate floating-point variables to the MC6888 l’s internal extended- 
precision registers, which are more precise than the single- or double-precision 
variables declared in FORTRAN, C, or Pascal. This means that, for instance, 

X = expression 
if (X .eq. 1.0) ... 

v J 
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may execute differently depending on whether X is allocated to an extended- 
precision register. If not, then the value of "expression" is rounded to the 
declared precision of X prior to the comparison, possibly changing its value; if X 
were so allocated, the rounding would not occur, affecting the subsequent com- 
parison. 

This optimization may adversely affect some programs that depend on the store 
to force rounding to storage precision. Those programs can be compiled without 
optimization, or with optimization and the additional option -f store, which 
insures that rounding to storage precision occurs when specified by a source 
language assignment 

-f store is implemented in Release 3.2 by not allocating floating-point vari- 
ables to extended-precision registers, -f store only has effect when 
-f 6 8 8 8 1 and -0 or -P are specified. 

The combination -f store -0 is safer than -0 in that a few programs either 
won’t work or work poorly with just -0. -0 is usually slightly more accurate; 
the difference in speed varies considerably among programs. 

The -f store option is not intended to affect anonymous temporary expres- 
sions. Thus the statement 



X = (x + y) - y 

V , 



may still be evaluated in extended precision until the final store to x, even if 
-f store is specified. This may cause problems with some artful programs, 
such as some in die Elef unt test suite of Cody and Waite. The Elefunt 
programs test elementary functions by computing various identities involving 
those functions; inaccuracies in the elementary function implementations are 
reflected by the extent to which the identities fail to be satisfied. To be meaning- 
ful, the arguments of the functions must be correctly related. That is, if x and 
x+k are both arguments to functions in the identity, then x must be such that 
x+k is computed exactly. Any error in x+k prior to evaluating the function 
under test will likely contaminate the results to such an extent that they are 
meaningless for measuring the errors in the function under test x is purified by 
perturbing it slightly until x+k can be computed exactly. 

A specific example in the Elefunt test of alog is an identity wherein x and 
y= ( 1 7 / 1 6 ) *x are both required. To compute y exactly on binary, decimal, 
and hexadecimal machines, over the limited range sqrt ( 1 / 2 ) <= x <= 
15/16, it suffices that x’s least significant four bits be zero; this is accom- 
plished by the following code: 

'l 

eight = 8.0 

x = (x + eight) - eight 
y=x+x/16.0 

s _J 



This code works correctly unless the expression (x + eight) - eight is 
evaluated in higher precision than the storage precision of x. If the latter occurs, 
then the rounding of (x + eight ) will not occur where intended and x will 
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not be purified. Since it is the intermediate expression (x + eight) that is 
of interest, -f store will not affect the result You must change the source 
code to 

'l 

eight = 8.0 

x = x + eight 
x = x - eight 
y = x + x / 16.0 

J 



-f single and 



-f single2 



Constant Expressions 



Then -f store will force rounding to storage precision after (x + eight), 
and y will be ultimately computed without any rounding error. 

Extended-precision expression evaluation is usually desirable because it minim- 
izes rounding error, and it is good practice to write out expressions with minimal 
stores of intermediate results. But programs which intend to force rounding 
errors to occur at certain times must sometimes be rewritten in this way when 
using extended precision. 

In C, floating-point expression evaluation and parameter passing are defined to 
occur in double precision, even if all operands are single-precision, cc’s 
-f single and -f single2 options permit more FORTRAN -like floating- 
point expression evaluation. 

-f single causes floating-point operations to be performed in single precision 
when operands are single-precision, -f single2 causes single-precision pro- 
cedure parameters to be passed as single-precision rather than double-precision. 
Do not mix C modules compiled with -f single2 with modules compiled 
without it. 



In particular, the following will not work if compiled with -f single 2: 



f — 


#include <math.h> 






float x,y; 
y = sqrt (x) ; 




V 




J 



Standard C libraries expect double-precision floating-point arguments. 

Constant expressions may be evaluated at compile time or at run time. From the 
point of view of the IEEE Standard, run-time evaluation of constant expressions 
is desirable in order to observe the rounding modes in effect at runtime and to 
preserve the IEEE exceptions. 

Sun’s FORTRAN compiler (f 77) doesn’t evaluate floating-point constant expres- 
sions at runtime, while cc and pc both do, using the default round to nearest 
mode. All three compilers evaluate integer constant expressions at runtime. 



Asun 

V microsystems 



Revision A of 12 September 1986 







Floating-Point Numerics 



Floating-Point Numerics 13 

2.1. Rounding Errors and Different Numerical Results 13 

Rounding Errors 13 

Different Numerical Results 13 

Detecting 111 Condition and Instability 14 

2.2. An Example- SPICE tdo 14 

2.3. Floating-Point Programs 15 

Assembly Language 15 

2.4. Precision in Decimal Digits 17 

Single Precision 17 

Double Precision 17 

Extended Precision 17 

2.5. FCVS Tests 17 

2.6. Floating-Point and Signal Handlers 19 

set jmp/long jmp 20 

Floating-Point Signals 20 

Signal #7 (SIGEMT) 20 

Signal #8 (SIGFPE) 20 

Why Not Signal by Default? 21 

2.7. Debuggers and Floating Point 21 

Signals 21 

Disassembly 22 

Printing all MC6888 1 Registers 23 



Printing all FPA Registers 23 

Printing in Single-Precision 23 

Printing in Double-Precision 23 

Printing in Extended-Precision 23 

Modifying Floating-Point Registers or Memory Locations 23 

2.8. FPA Recomputation 24 

Rounding Modes 25 

Exceptions 25 

Inexact Exceptions 25 

Traps 26 

Performance 26 



Floating-Point Numerics 



2.1. Rounding Errors and 
Different Numerical 
Results 



Rounding Errors 



Expression evaluation and precision are important because floating-point calcula- 
tions are usually inexact. The set of numbers representable in a particular finite- 
precision format is not closed under the usual arithmetic operations. When the 
computed result differs from the exact infinite-precision result, we say that 
roundoff or a rounding error has occurred. Correctly implemented floating point 
delivers the exact result if it is representable, or one of at most two nearest 
representable numbers otherwise. Computers vary as to the mle used for choos- 
ing between the two nearest neighbors; some computers with defective arithmetic 
do not guarantee to produce representable results correctly or to produce one of 
the nearest neighbors for unrepresentable results. 

Sun’s computers generally follow the IEEE Standard for Binary Floating-Point 
Arithmetic, which specifies that, by default, inexact results return the closer of 
the two nearest neighbors, and if they are equally close, the one with the least 
significant bit 0. This mle is the default "round to nearest" rounding direction 
mode. Rounding errors caused by the correct application of the rounding mle to 
inexact results are inevitable consequences of finite-precision arithmetic, rather 
than defects in the machine. 



Different Numerical Results 



Different floating-point implementations produce different results. Even when all 
floating-point operations conform to the IEEE Standard, elementary transcenden- 
tal functions will usually vary depending on the hardware. Unlike the rational 
operations (+, -, *,and /) and the algebraic function sqrt, it is not econom- 
ical to produce correctly rounded results for the elementary exponential, loga- 
rithmic, and trigonometric transcendental functions. Consequently different 
implementations make different trade-offs between speed and accuracy. 

Differences arise even if no transcendental functions are involved. The 
MC6888 1 and Sun FPA both conform to the IEEE Standard, yet they produce 
different results because the MC6888 1 uses extended- precision registers, which 
are not present on the FPA. The MC6888 l’s results are usually more accurate. 
The MC6888 l’s results can vary depending on which of the -O, -P,or 
-f store flags are used, since these affect how often results are rounded from 
extended precision to storage precision. Finally, results computed with the A79J 
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Detecting 111 Condition and 
Instability 



2.2. An Example- SPICE 
tdo 



mask set MC6888 1 may differ from those computed with the later versions of the 
MC6888 1. 

Sun’s software floating-point conforms to the IEEE Standard’s default modes 
and agrees with the Sun FPA on operations specified by the Standard in the 
default modes. 

Due to performance requirements and the fact that its microcode space is com- 
pletely full, the Sky FFP does not always produce the results specified in the 
IEEE Standard, even for the rational operations +, *, and /. For these 

operations and sqrt, computed results sometimes differ by one bit from the 
standard results; at other times the signs of zero results vary from the Standard. 

If you specify cc’s -f single or -fsingle2 options, results may be dif- 
ferent from what they would have been otherwise. 

When the answer to a problem is inordinately sensitive to changes in the input 
data, it is said to be ill-conditioned. Ill condition is a property of a problem and 
input data and is not affected by changing the algorithm used to obtain the result. 

Even for well-conditioned problems, some numerical algorithms are much more 
sensitive to roundoff than others. The sensitive algorithms are said to be 
unstable. Instability can sometimes be cured by changing algorithms. 

The distinction between ill-condition and instability is often academic in large 
realistic applications too complicated to be analyzed. The symptoms are the 
same: a computation previously believed to be correct on another machine, or 
with a different compiler or operating system, or even with different compiler 
options, now gives drastically different results with the same input data. Natur- 
ally, such results suggest that the new hardware, compiler, or operating system is 
defective. How can that possibility be separated from the effects of ill condition 
or instability? 

One way to investigate sensitive results is to introduce random small changes in 
the input data and observe the effect on the results. But the right random changes 
may be difficult to determine or awkward to insert. The four rounding direction 
modes specified in the IEEE Standard allow one to obtain similar effects. 

On Sun-3 ’s with an MC6888 1 installed, you can vary the IEEE rounding modes 
with the assembly-language function, _f pmode_, mentioned in "How to Use 
Special Features of the IEEE Standard" in Chapter 3. Run a program first in the 
default round-to-nearest mode, then in another mode such as round-toward-zero. 
Compare the results, and if numerical results vary enough to bother you, then you 
are either trying to solve an ill-conditioned problem or using an unstable algo- 
rithm. 

SPICE is an integrated circuit emulation program. One of the SPICE 2G.6 
input files listed in Appendix E, tdo, models a tunnel diode oscillator. In the 
transient analysis section of the program’s output, a voltage vl is monitored as 
a function of time. The computed behavior is that the voltage is initially a con- 
stant 0. 12 v, then either declines abruptly to about 0.02 v or else jumps equally 
abruptly to about 0.48 v. Which behavior is computed seems to depend on ran- 
dom factors of code generation causing slight variations in rounding errors. For 
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2.3. Floating-Point 
Programs 



instance, the abrupt decline is obtained with -fsoft -0 and -ffpa -O, 
while the abrupt rise is obtained with -fsky -0 and -f68881 -0,inall 
cases compiling with software Release 3. 1 and default roundings. Experiments 
with -f68881 -0 show that changing the rounding direction mode from 
round to nearest to round toward zero, or changing the rounding precision mode 
from extended to double, is sufficient to change the computed behavior. Since 
oscillators usually exhibit positive-feedback physical characteristics, it seems 
likely that computing the transient behavior is an inherently ill-conditioned prob- 
lem. On the other hand, determining the steady-state oscillation frequency might 
be well-conditioned. 

Sun recommends that you write programs using floating-point in higher-level 
languages whenever possible. 

FORTRAN is the best-suited compiled language that Sun supports for those parts 
of large applications where floating-point usage is substantial, and data struc- 
tures, control structures, and input and output are simple. FORTRAN has more 
intrinsic types, operators, and functions designed to support floating-point com- 
putations and ease writing correct, compact, easily optimized numerical pro- 
cedures. 

Many other desirable language features are poorly designed in FORTRAN or 
lacking altogether, so Sun makes it relatively easy to call FORTRAN computa- 
tional modules from programs written in C or Pascal. 



In Pascal, you just declare such modules to be "external fortran". 
In C, remember that a FORTRAN declaration like 




usually corresponds to a C declaration like 




Assembly Language 



For more about programming on Sun’s computers in FORTRAN and Pascal, see 
the Sun FORTRAN Programmer’ s Guide and the Sun Pascal Programmer’ s 
Guide. 

The assembly language interfaces to -f switch, -fsoft, and -fsky float- 
ing point are not documented and vary in different software releases. If neces- 
sary, you can examine compiled code generated by the -S option to infer the 
interface conventions in effect for a particular release. 

For -f 6 8 8 8 1 and -ffpa, the interface is more permanent since specific 
features have been added to the Sun-3 assembler to accommodate these floating- 
point units. The following lines of code are comparable: 



A 
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Source: 

— 
x(i) = x (i) + c * y(i) 

/ 



MC6888 1 assembly code: 



r 






s 




fmoved 


c,fpO 






fmuld 


ay6+, fpO 






faddd 


ax6, fpO 






fmoved 


fpO, ax0+ 




V 






j 



FPA code: 

. 

fpmoved c,fpaO 
fpmuld ay@+,fpaO 
fpaddd ax@,fpaO 
fpmoved fpaO,ax0+ 

s > 



In both examples, ax and ay are pointers to x(i)and y (i) , respectively. 

The syntax of MC6888 1 instructions generally follows that of Motorola’s 
MC6888 1 manual, with the following exceptions: 

□ Instructions are lower case 



□ No dot separates the opcode from the operand type 

□ CPU address mode notations are different 

Furthermore, unlike examples in the MC6888 1 manual, dyadic operations 




FPA instruction syntax is designed to parallel that of MC6888 1 instructions 
where it makes sense to do so. Op codes begin with "f p" rather than "f ' " and 
register names are fpaO . . fpa31 rather than fpO . . fp7. Register-to- 
register operations have type s or d rather than x. 
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2.4. Precision in Decimal 
Digits 



Single Precision 



Double Precision 



Extended Precision 



2.5. FCVS Tests 



The relationship between binary floating point and decimal floating point is one 
of the most widely misunderstood aspects of floating-point arithmetic. No binary 
floating-point format has an exactly equivalent number of decimal digits, since 
the relative spacing of representable floating-point numbers varies in different 
places for different bases. The following describe the IEEE formats: 

IEEE single precision is always more precise than 6 decimal digits and always 
less precise than 9 decimal digits. Seven decimal digits are often less precise but 
sometimes more; eight decimal digits are often more precise but sometimes less. 

IEEE double precision is always more precise than 15 decimal digits and always 
less precise than 17 decimal digits. Sixteen decimal digits are sometimes less 
precise and sometimes more. 

MC6888 1 extended precision is always more precise than 18 decimal digits and 
always less precise than 21 decimal digits. Nineteen decimal digits are often less 
precise but sometimes more; 20 decimal digits are often more precise but some- 
times less. 

The following examples illustrate the confusion surrounding precision issues. 
The Federal Compiler Validation Service (FCVS) has a series of programs 
intended to measure conformance of implementations to the FORTRAN standard, 
ANSI X3.9-1978. While investigating why certain test programs failed on Suns 
in mysterious ways that depended on the setting of the -f . . . and -O options, 
it became apparent that the programs were setting standards that could only be 
passed by luck. Here are two examples taken from the FCVS78 Version 2.0; 
official FCVS versions since 1983 have been modified to display the message 
"INFORMATIVE - UNDER REVIEW BY FSTC" for certain tests in programs 
371..374, 818, and 820. 

The following fragment is from test 10 of FCVS program 374, intended to test 
single-precision tangent: 




This fragment is intended to be a test of tan, to see if 




which is approximately 1024 - .0003. So the passing criterion seems rather gen- 
erous: 




But whether this test passes or fails is a matter of random luck for single preci- 
sion. Suppose that the only error is the error in the value assigned to P IVS; it 
was intended to be Jt, but since K is not representable, the actual value contained 
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in PI VS is not n but Jt+e. Conventional error analysis shows that the passing 
criterion is equivalent to: 

< 

abs(e) < (2/3) * (0.1) * 2**-20 = 6.4e-8. 

V > 



But binary floating-point arithmetic of p significant bits cannot in general 
approximate a number to better than half a unit in the last place, or 2**(1 -p) for 
numbers between 2 and 4; IEEE single precision has 24 significant bits, so arbi- 
trary numbers between 2 and 4 would not be expected to be approximated to 
better than e = 1.2e-7, and in particular, the IEEE single-precision number closest 
to 7t differs from 7t by £ = 8.7e-8. Thus there is no hope of passing unless addi- 
tional rounding errors are made in the multiplication or tangent; these additional 
rounding errors must be aptly chosen to cause the final result to be acceptable! 



The following fragment is from test 12 of FCVS program 820, intended to test 
complex cosine: 



COMPLEX AVC, BVC 
REAL R2E (2) 

EQUIVALENCE (AVC, R2E) 





AVC 


= CCOS ( ( 




IF 


(R2E (1) - 


40121 


IF 


(R2E ( 1 ) - 


40122 


IF 


( R2 E ( 2 ) + 


40120 


IF 


(R2E (2 ) - 



3.1416, 0.0) 
0 . 99725E+00) 
0 . 99736E+00) 
0 . 50000E-04) 
0 . 50000E-04) 



* (- 10000 . 0 , 0 . 0 )) 
fail, 40122, 40121 
40122, 40122, fail 
fail, pass, 40120 
pass, pass, fail 



This is intended to test whether 

> 

ABS (COS ( le4 * 3 . 1416) -. 997305) < 5 . 5e-5 
v ; 



The correct result cos(3 1416.0) is .9973027. But once again, passing is a matter 
of random luck on binary computers. Suppose the only error is the conversion of 
3. 1416 to a machine-representable value 3. 1416+e. Error analysis shows that the 
passing criterion is equivalent to 

> 

-8 . Oe-8 < £ < 7 . Oe-8 

^ / 

but as before, the best general bound for the approximation error is 1.2e-7, and in 
particular the closest IEEE single-precision number to 3.1416 differs from it by 
l.le-7. Machines with 24 significant bits or less can only pass the indicated cri- 
terion by lucky rounding in the multiplication and cosine. 

Oddly enough, this same test program considers acceptable a nonzero imaginary 
part of a complex cosine of a real argument. 
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2.6. Floating-Point and 
Signal Handlers 



UNIX signals are asynchronous events, and when a signal handler uses the same 
floating-point device as the process it interrupted, confusion can ensue. Conse- 
quently part of Sun’s signal handling protocol is designed to save and restore the 
state of the floating-point device currently in use, but it is good practice to design 
signal handling routines to be as short as possible and to avoid inessential 
floating-point computations. 

Signals, once generated, are passed by UNIX to the interrupted process through a 
routine named _s igt r amp which calls a signal handler defined by the user 
with sigvec(2), signal (2), or signal (3f). 

Before calling the user handler, _s igt ramp saves certain CPU and floating- 
point state on the stack; if the signal handler terminates normally by returning to 
_s igt ramp, the CPU and floating-point state are restored from the saved state. 
The CPU state includes the contents of the program counter (pc), stack pointer, 
status register, signal stack and mask information, and registers aO-al and 
dO-dl. The floating-point state depends on which signal occurred and what kind 
of floating-point hardware is in use. In the case of signal #8, SIGFPE, no 
floating-point state is saved or restored; for all other signals, the hardware state of 
the Sky FFP, MC6888 1, or Sun FPA is saved when the signal occurs and restored 
when the signal handler returns. 

Note that just as the signal handler is expected to follow the C compiler’s con- 
ventions and restore CPU registers a2-a7 and d2-d7 if it uses them, it is also 
expected to restore certain floating-point registers, as described below. 

In the case of the Sky FFP, the complete state is saved, consisting of the contents 
of the program counter of the Sky FFP and the contents of the four double- 
precision floating-point data registers. Before starting the signal handler, the Sky 
FFP is reset to accept new instructions. The contents of the data registers are 
unchanged. On return from the signal handler, the Sky FFP program counter and 
data registers are restored to their state prior to interruption. 

In the case of the MC68881, the registers saved and restored include fpO-fpl 
and fpcr/ fpsr/fpiar. The signal handler is responsible for saving and res- 
toring fp2-fp7 if ituses them. 

In the case of the Sun FPA, the registers saved and restored include the 
MC68881’s fpO-fpl and fpcr/fpsr/fpiar and the FPA’s instruction 
pipe, fpastatus, fpamode, imask, load_ptr, ierr, and fpaO- 
fpa3. Before calling the signal handler, imask is set to 0, fpamode 
is set to 2 (the default), and f psr and fpcr are set to 0. The signal handler is 
responsible for saving and restoring fpa4-fpa31 if used. Upon returning 
from the signal handler, all the control and data registers are restored to their 
values at the time of the signal. 

The exact format of the saved floating-point state is subject to change and should 
not be relied upon. Programs that use nonstandard returns from signal handlers 
may not work under new software releases. 
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set jmp/ long jmp 



Floating-Point Signals 



The code implementing set jmp and long jmp does not save or restore any 
floating-point state. Use great care in creating programs that involve these func- 
tions and floating-point arithmetic. 

Two floating-point signals, SIGEMT and SIGFPE, are used by software to test 
for the presence of an MC6888 1 and as signals for some error or exception con- 
ditions. 



Signal #7 (SIGEMT) 



Programs compiled with -fswitch, -f68881,and -ff pa test for the pres- 
ence of an MC68881 by performing an MC68881 instruction that initializes the 
f per and fpsr registers to their proper IEEE defaults. If an MC6888 1 is 
present, the instruction executes uneventfully. If no MC6888 1 is present, a 
SIGEMT signal is generated and handled by a signal handler installed for the 
duration of this test. The SIGEMT signal also affects use of the debuggers 
adb, dbx and dbxtool; see "Debuggers and Floating Point" below. 



Signal #8 (SIGFPE) 



SIGFPE signals certain conditions that might be errors, exceptions, or opportun- 
ities. Generally these conditions arise synchronously during arithmetic opera- 
tions in the underlying process, but they are infrequent enough that they are most 
efficiently treated as if they were asynchronous events. Since they usually tell 
something about the state of the floating-point computation, the floating-point 
status is not saved, reset, and restored as it is for all other signals. 



Not all SIGFPE signals have to do with floating point. This is a complete list 
of defined SIGFPE codes: 



# define 


FPE_INTD I V_TRAP 


0x14 


/* 


#def ine 


FPE__CHKINST_TRAP 


0x18 


/* 


#def ine 


F P E_T RAP V_T RAP 


Oxlc 


/* 


#def ine 


FPE_FLTBSUN_TRAP 


OxcO 


/* 


#def ine 


FP E_F LT I NE X_T RAP 


0xc4 


/* 


#def ine 


FPE__FLTDIV__TRAP 


0xc8 


/* 


#def ine 


FP E_F LTUND_TRAP 


Oxcc 


/* 


#def ine 


FPE__FLTOPERR_TRAP 


OxdO 


/* 


#def ine 


FP E_FLTOVF_TRAP 


0xd4 


/* 


#def ine 


FPE_FLTNAN_TRAP 


0xd8 


/* 


#def ine 


FPE__FPA_ENABLE 


0x400 


/* 


#def ine 


FP E_F P A_E RROR 


0x404 


/* 



integer divide by zero */ 

CHK [CHK2] instruction */ 

TRAPV [cpTRAPcc TRAPcc] instr */ 
[floating branch | set on unordered ] */ 
[floating inexact result] */ 

[floating divide by zero] */ 

[floating underflow] */ 

[floating operand error] */ 

[floating overflow] */ 

[floating Not-A-Number] */ 

[FPA enable reg not set] */ 

[FPA bus error] */ 



Of these, INTDIV, CHKINST, and TRAPV are related to integer arithmetic or 
non-arithmetic instructions. FLTBSUN . . . FLTNAN are generated by the 
MC6888 1 and can be indirectly generated by the FPA. FPA_ENABLE and 
FPA_ERROR can only be generated by the FPA. 

The FP E_FP A_ERROR code is used by the FPA to indicate that it cannot per- 
form a particular operation; this may be because the operation is not a valid FPA 
instruction or the data is outside the domain accepted by the Weitek chips or the 
FPA microcode. The latter exceptions are part of normal FPA operation. 

SIGFPE must be enabled and assigned to a signal handler that recognizes the 
FPE_FPA_ERROR code and handles it properly when an FPA is to be used. 
(This is done automatically by -ffpa). If SIGFPE is subsequently totally 
disabled, an FPE_FPA_ERROR signal will cause an infinite loop because the 
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Why Not Signal by Default? 



2.7. Debuggers and 
Floating Point 

Signals 



FPA is never reset and cannot progress past the troublesome instruction. 

Floating-point S IGFPE codes do not arise for programs compiled with 
-f sof t or -f sky. They do not arise for programs compiled with -f 68 8 8 1 
unless IEEE traps are enabled on the MC6888 1 by setting bits in the Exception 
Enabled byte of the fpcr register. The MC6888 1 exception enabled bits also 
affect FPA results during recomputation; see "FPA Recomputation" below. 

Users new to the IEEE Standard often wonder why seemingly catastrophic events 
such as division by zero or overflow do not cause immediate error termination of 
the program. They have been conditioned to expect the latter by various systems 
of arithmetic that were unable to provide any reasonable response for exceptional 
conditions. 

The IEEE Standard, however, provides default responses for floating-point 
operations that are meaningful and adequate for many purposes: Infinities and 
NaNs are the usual results of overflow and division by zero. Sun’s floating point 
meets the standard by silently providing the correct default result. This is usually 
satisfactory when the result of an operation is a floating-point quantity; other 
cases, sometimes troublesome, include comparisons, conversion to integer for- 
mat, and conversion from binary to ASCII. Programmers who haven’t con- 
sidered the possibility of floating-point operands that might be infinities or NaNs 
sometimes get surprising results in these cases. 

When an MC6888 1 is available, the S IGFPE signals FPE_FLTBSUN_TRAP 
through FPE_FLTNAN_TRAP can be enabled to occur when specified floating- 
point exceptions arise. For instance, to cause FPE_FLTOPERR_TRAP to be 
signalled whenever an invalid operation occurs such as converting a NaN, 
infinity, or too-large value to an integer, it is necessary to set the OPERR bit in 
the Exception Enable Byte in the MC68881’s fpcr register, described in 
Motorola’s MC6888 1 User’s Manual. This can be done either in assembly 
language, or by using the _f pmode_ subroutine described in the section "How 
to Use Special Features of the IEEE Standard" in Chapter 3. 

Floating-point hardware or software use will affect the way that you work with 
debuggers on Sun systems. Signal handling and disassembly operations are par- 
ticularly affected. For more information about debuggers and floating-point, see 
Appendices B and C in this manual. 

When using the debuggers adb, dbx, and dbxtool, be aware that certain 
signals occur in the normal course of floating-point operations. As mentioned 
above, a MC6888 1 instruction is attempted to determine if an MC6888 1 is 
present If no MC6888 1 is present, S I GEMT is signalled and is caught by the 
debuggers, stopping the process being debugged. You can then either continue 
or, if you prefer, instruct the debugger in advance to ignore the signal by issuing 
the command 7 : i in adb or ignore 7 in dbx. This command passes the 
signal directly to the debugged process, bypassing the debugger. 

When an FPA is installed, S IGFPE is signalled each time FPA recomputation is 
required, as discussed below. When these signals occur, you may either con- 
tinue, or issue 8 : i or ignore 8 to ignore them. Almost every process that 
uses the FPA requires recomputation for at least the first inexact result. 
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Disassembly 



MC6888 1 instructions can always be disassembled using MC6888 1 assembly 
language mnemonics. FPA instructions, however, are simply move 1 instruc- 
tions to particular addresses, and consequently cannot be positively distinguished 
from other move 1 instructions. To allow control over this disassembly, adb 
and dbx provide ways to specify whether or not movel instructions should be 
considered as potential FPA instructions. 



Sun’s compilers generate FPA instructions using movel’ s with "long absolute" 
and "address register with displacement" effective address modes. For example: 





fpadds 


<eadatal>, fpaO 


N 




fpmoves@l 


f paO , <eadata2> 




V 






J 



generates exactly the same code as 



r 


movel 


<eadatal>, 0xe0000380 


> 




movel 


al@ (OxcOO) , <eadata2> 




V 









The debuggers can recognize two patterns of move instructions as FPA instruc- 
tions. First, instructions that access the FPA using absolute addressing may be 
recognized simply by examining the address; second, instructions that access the 
FPA using base + 16-bit displacement addressing may be recognized if the base 
register used has been designated as a pointer to the FPA’s virtual address. 

Note that you can make this designation at runtime, independent of the actual 
register contents. The default is to recognize absolute-mode FPA instructions 
only if an FPA is present. Base-mode FPA instructions are assumed to be moves 
by default, since no register is assumed to be an FPA base register. 

FPA disassembly is enabled in adb by the command 1>F and disabled with 
0>F. 1>B designates al as an FPA base register, while -1>B indicates that 

no register is being used as an FPA base register. The first example above would 
be disassembled as an FPA instruction if F were 1; the second if F were 1 and 
B were 1. Otherwise they would be disassembled as movel’s. 

In dbx, 
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dbxenv fpabase off 



Printing all MC68881 
Registers 



Printing all FPA Registers 



Printing in Single-Precision 



Printing in Double-Precision 



Printing in Extended- 
Precision 



Modifying Floating-Point 
Registers or Memory 
Locations 



indicates that no register is assumed to be an FPA base register. 

In adb, $R displays MC6888 1 registers fp0-fp7 in hexadecimal (hex) and 
extended-precision decimal, and fpcr, fpsr, and fpiarinhex. 

In dbx, the command & $ f p 0 / nE displays n MC6888 1 registers starting with 
fpO, in hex and extended-precision decimal. 

In adb, $x displays FPA registers fpaO-f pal 5 in hexadecimal, single- 
precision decimal, and double-precision decimal. It displays fpamode, 
fpastatus, state, imask, ldptr, ierr, and instruction pipe registers 
in hex. $X displays FPA registers fpal6-fpa31 in hex, single decimal, and 
double decimal, and fpamode, fpastatus, state, imask, ldptr, 
ierr, and instruction pipe registers in hex. 

In dbx, the command &$fpaO/nF displays n FPA registers starting with 
fpaO, in hex and double-precision decimal. To display FPA registers in hex and 
single-precision decimal, use &$fpaOs/nf . 

In adb, a 4-byte memory location, f pan, or a register dn can be displayed in 
single-precision decimal; use addr/f or <fpan=f or <dn=f in adb. 

In dbx, use addr / f to print a single location in hexadecimal and single- 
precision decimal, addr/nf to print n locations. The notation &$regname 
can be used as a handle for printing registers; note that registers do not actually 
have addresses. 

In adb, an 8-byte memory location, f pan, or a register pair dn : dn+1 can be 
displayed in double-precision decimal; use addr/F or <fpan=F or <dn=F. 

In dbx, use addr /F to print a memory location or register in hexadecimal and 
double-precision decimal, addr /nF to print n registers or memory locations. 

A 12-byte memory location or fpn can be displayed in extended-precision 
decimal: In adb, use addr/eor <fpn=e. 

In dbx, use addr/Eor &$fpn/E. 

In adb you cannot load memory or registers with data in an ASCII floating- 
point representation. Hexadecimal data can be used to load one 32-bit word at a 
time. The words of a register can be designated separately, to make loading 
MC6888 1 and FPA registers easier. Thus the most significant, middle, and least 
significant words of f pi can be referred to as fpl, fql, and frl. The 
more significant and less significant words of the FPA register f pal can be 
referred to as fpal and fqal. 

In dbx you can set a register or memory location with a value expressed in 
floating-point decimal notation. 

dbx converts data into extended-precision format for an MC6888 1 register, and 
into double- or single-precision format for an FPA register, depending on 
whether a name in $fpaO. .$fpaO or $fpaOs. .$fpa31s is used. The 
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data is converted according to the declared type of the variable for variables 
defined in a compiled language, dbx also provides the command 




to set an MC6888 1 register to an arbitrary 3-word bit pattern, formed by con- 
catenating the 32-bit values x, y and z. No such command exists for FPA regis- 
ters in Release 3.2. 



2.8. FPA Recomputation 



The Sun FPA relies on the Sun-3 CPU’s MC6888 1 to perform many operations 
that would be difficult or inefficient to perform with the Weitek 1 164/1 165 chip 
set. These operations include any in which IEEE arithmetic exceptions arise. 
When one of these exceptions is detected in the Weitek chips, the operation is 
terminated and a bus error notifies the CPU that recomputation on the MC6888 1 
is required. 

The operating system intercepts the bus error, diagnoses it as coming from the 
FPA, and sends a SIGFPE to the process in question. That process is responsi- 
ble for having a correct S IGFPE signal handler installed that determines if the 
SIGFPE is due to an FPA request for recomputation and handles it accordingly. 
When -f fpa is requested, such a signal handler is installed by default; the fol- 
lowing is a simplified version of that default handler. 



♦include <signal.h> 



sigfpe_handler (sig, code, scp) 
int sig, code; 

struct sigcontext *scp; 



if (code == FPE_FPA_ERROR) { 
if (fpa_handler (scp) ) 
return; 

1 

} else { 

fprintf (stderr, "SIGFPE Handler pc %X code %X \n", scp->sc_pc, code); 
f flush (stderr) ; 



abort ( ) ; 



fpa_handler ( ) is the standard FPA recomputation routine contained in 
libc. It returns 1 if the FPA error was fixed by recomputation and 0 otherwise, 
as in the cases that the FPA error was caused by hardware failure or incorrect 
software. fpa_handler ( ) is used to support the IEEE Standard arithmetic 
features pertaining to rounding modes, exceptions, and traps. 
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Rounding Modes 



Exceptions 



Inexact Exceptions 



fpa_handler ( ) uses the MC6888 1 to recompute the desired result. The 
MC6888 l’s single and double rounding precision modes are used to obtain the 
results anticipated by the IEEE Standard for simple single- and double-precision 
operations. More complex instructions, such as elementary transcendental func- 
tions and combination instructions, such as add followed by multiply or multiply 
followed by add, are recomputed in extended rounding precision mode to obtain 
maximum accuracy. 

Rounding modes are supported directly by the Weitek 1164/1165 for simple 
instructions; the elementary transcendental functions, however, are microcoded. 
When an elementary transcendental function is called with a rounding mode 
other than the default (round to nearest) recomputation is invoked in order to 
compute the function on the MC68881. 

Exceptions are reported by the Weitek chips after each operation, but those chips 
do not accumulate exceptions or allow optional trapping on them. Consequently, 
whenever the Weitek chips report one of the exceptions other than inexact, the 
FPA hardware signals that recomputation is required. The CPU receives that sig- 
nal when it attempts to perform a subsequent FPA instruction; the CPU’s pro- 
gram counter points to the subsequent instruction rather than to the instruction 
which suffered the exception; the latter’s address is lost. When an exception is 
signalled, the result delivered by the Weitek chips is not written to the destina- 
tion FPA register specified by the FPA instruction, since in many cases that 
would destroy input data required to recompute the correct result 

The recomputation routines examine the FPA’s instruction pipe, clear it, and 
recompute the first instruction in the pipe, which is the one that generated die 
exception. The recomputed result is placed in the FPA register that was the 
instruction’s intended destination, and any IEEE exception flags raised during the 
course of the computation are OR’ed into the MC6888 l’s Accrued Exception 
byte in the f psr register. 

After recomputation occurs, the contents of the FPA’s f pastatus register and 
the MC6888 l’s Condition Code, Quotient, and Exception Status bytes are 
undefined, except that after recomputation of an FPA comparison instruction 
(which happens when one or both operands are NaNs), f past atus is properly 
defined so it can be read and then written into the MC68020’s status register. 

Inexact exceptions arise more often than not during floating-point operations. If 
inexact is the only exception that arose, the Weitek-computed numerical result is 
correct. Causing a recomputation on each such exception would unacceptably 
degrade performance. Consequendy the FPA has a register called imaskthat 
allows software to control whether inexact exceptions invoke recomputation. 
Sun-supplied software uses imask in conjunction with the MC68881 ’s status 
and control registers to maintain an accurate record of inexact exception 
occurrences without recomputing for each one. Sun’s software conventions for 
imask use may be defeated if users alter imask other than indirectly through 
the function _f pmode_ described below. 

When a -f fpa program begins execution, imask is set to 1 so that the FIRST 
inexact exception detected by the FPA causes a recomputation. After the recom- 
putation, the inexact Accrued Exception bit in the MC68881’s fpsr is turned 



m sun 

XT microsystems 



Revision A of 12 September 1986 




26 Floating-Point Programmer’s Guide 



Traps 



Performance 



on and imask is cleared UNLESS the inexact operation Exception Enable bit is 
set in the MC6888 l’s f per . The latter bit, normally cleared, may be set by the 
user to indicate that EVERY inexact FPA operation should cause recomputation 
by the MC6888 1, with a trap occurring if the MC6888 1 recomputation is also 
inexact 

If you use either the function _f pstatus_ to clear the MC68881 fpsr’s 
inexact Exception Accrued bit, or _f pmode_ to set the fper’s inexact opera- 
tion Exception Enabled bit, then the FPA’s imask will be set so that the next 
inexact operation will cause recomputation. 

Traps on FPA instructions can be obtained by setting bits in the MC6888 l’s 
Exception Enable byte in the f per register. When an exception is detected by 
a Weitek chip the FPA calls for recomputation; if that recomputation causes a 
MC68881 exception corresponding to a set Exception Enable bit then a 

MC68881 trap occurs, manifested as a SIGFPE with a code equal to one of 
— 

FPE_FLT { BSUN, INEX, DIV, UND , OPERR, OVF , NAN }_TRAP 




The MC68881 trap occurs after the FPA recomputation takes place; the pc in the 
SIGFPE handler parameter sep will point to an FPA instruction in the user’s 
code, the same one pointed to by the orgignal recomputation request The 
address of the last MC6888 1 instruction executed in fpa_handler ( ) is in 
theMC68881’s fpiar register. 

Writing a SIGFPE trap handler for an MC6888 1 trap generated indirectly by 
the FPA requires considerable care. The same trap handler must handle both the 
FPA’s request for recomputation and the MC6888 l’s trap signal. 

Sun does not support user-written replacements for f pa_handler ( ) . If you 
choose to write one, be sure to meet all hardware error-handling requirements; in 
particular, clear the FPA’s pipe before attempting any other FPA instructions. 
Otherwise an infinite loop of SIGFPE’s is likely to result. 

Because the treatment of floating-point exceptions is so variable, most floating- 
point programs that run reliably on machines with different types of arithmetic 
do not cause any floating-point exceptions to occur other than inexact. Conse- 
quently Sun’s FPA recomputation scheme was designed without consideration of 
efficiency other than in the inexact case. 

Programs that generate many exceptions may run more slowly than software 
floating point due to the intervention of the operating system to handle bus errors 
and signalling. The most likely cases where this might arise are computations 
which frequently underflow or which handle many NaNs. If recomputation per- 
formance is an issue, it may be partially ameliorated by using the Weitek chips in 
"fast" mode. "Fast" mode, although inconsistent with the IEEE Standard, causes 
no ill effects on programs that work correctly on a variety of non-IEEE 
machines, since those machines probably underflow abruptly to zero without 
creating subnormal (denormalized) numbers as IEEE machines do. In the 
Weitek chips’ normal "IEEE" mode, recomputation is invoked if either operand 
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or the result is a subnormal number. In the Weitek chips’ "fast" mode, recompu- 
tation is only invoked if the result would be a subnormal number; subnormal 
operands are treated as zero without generating any exception. 



The Weitek f pamode register is intended for use by Sun software. If you alter 
it directly using an assembly language instmction 



fpmove #0x3,fpamode | Set "FAST" mode+standard rounding. 



or by an equivalent FORTRAN call 



< 




\ 




call fpamode(3) 




s 







be careful to preserve the values of mode bits that you do not intend to change. 
The default value of f pamode is 2, corresponding to rounding to nearest mode, 
treating subnormal operands in "IEEE" mode, and converting floating to integer 
in round-toward-zero mode. 



The following program illustrates performance aspects in the face of underflows 
— a few (c = 0.49) or a lot (c = 0.51). 




r 






A 


O 

II 

o 


O 

II 

o 

cn 




Options 


112 


491 


-0 


-f sof t 


100 


162 


-0 


-f68881 -fstore 


20 


70 


-0 


-f 68881 


15 


13800 


-0 


-f fpa 


15 

V 


15 


-0 


-ffpa (call fpamode(3) uncommented) 

j 
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IEEE-Standard Conformance 



The conformance of Sun software to the IEEE Standard varies depending on 
what compiler code generation options are chosen. The description in this 
chapter refers to the underlying arithmetic of Release 3.2 FORTRAN, Pascal, and 
C. Some of the features of the IEEE Standard that are unsupported in higher- 
level languages must be accessed through constructed subroutines, generally in 
assembly language; examples of these are listed later. 



3.1. Supported IEEE 
Standard Features 

Numeric Formats 



ASCII-to-Binary Conversion 



Rounding Direction Modes, 
Exceptions, and Traps 



Numerical Results 



All Sun floating-point options support single and double precision. Extended 
precision is supported only in assembly language for the MC6888 1, as are round- 
ing precision modes. Signalling NaNs are only recognized with -f 68 881 and 
-f fpa. 

Sun’s current software does not support IEEE-standard conversion from ASCII 
strings representing decimal numbers to binary floating point or the reverse. Such 
conversion is always performed in the default rounding mode, signals no excep- 
tions, and errors may be slightly larger than prescribed by the Standard, particu- 
larly in double precision. Assembly language programmers can obtain standard 
conversion using an A93N MC6888 1, if one is present. 

Rounding mode, exception, and trap support are handled differently by the vari- 
ous floating-point options: 

□ -f sof t and -f sky don’t support rounding modes, exceptions, or traps 

o -f68881 supports rounding direction and precision modes, exceptions, and 
traps 

□ -f fpa supports rounding direction modes and exceptions but not the full 
trapping mechanism recommended by the IEEE Standard 

Correct single- and double-precision numerical results are obtained for all opera- 
tions (except conversion between ASCII and binary) in the following cases: 

□ -fsoftin round-to-nearest mode 

□ -f fpa, provided the MC68881’s version is A93N 

□ -f 68881, if the MC68881 is A93N and MC68881 rounding precision is set 
to double 
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□ -f sky’s nonzero results may vary from the Standard’s prescription by one 
bit; zero results may have the incorrect sign. 

Correct results for -f f pa and -f 6 8 8 8 1 depend on the presence of an A93N 
mask version MC6888 1. 



MC68881 Rounding Precision 
Mode 



Extended-precision format is optional in the IEEE Standard, but if provided 
rounding precision modes must also be supplied. These modes govern the round- 
ing of extended-precision results; mode settings make it possible to round 
extended-precision results to shorter precisions in order to simulate results 
obtainable on systems without extended precision. 

The default MC6888 1 rounding precision mode is extended. This means that 
extended-precision operations upon extended-precision operands produce 
extended-precision results and exceptions that conform completely to the Stan- 
dard. The design of the MC6888 1 makes this mode substantially more efficient 
than single or double rounding precision modes. - f 6 8 8 8 1 sets up the default 

MC6888 1 settings, including extended rounding precision mode. 



However, a double-precision operation that computes a double-precision result 
from double-precision operands will then be computed on the MC6888 1 with 
two roundings: 



r 








fmoved 


x, fpO 


| This operation is exact. 




fmuld 


y, fpO 


| Round extended-precision 


fpO here. 


fmoved 


fpO, z 


| Round double-precision z 


here . 

j 



Ordered Comparisons 



3.2. How to Use Special 
Features of the IEEE 
Standard 



From the point of view of verifying double-precision conformance to the IEEE 
Standard, multiple roundings on a single arithmetic operation cause incorrect 
results to be generated occasionally, even though default extended-precision 
accumulation is generally faster and more accurate. The double rounding preci- 
sion mode can be used to insure that only one rounding occurs (in f mu Id in the 
example above), mainly for the purpose of verifying conformance with the IEEE 
test vectors and the Paranoia program discussed later. 

Programs should not generally use double rounding precision mode unless they 
require, on an operation-by-operation basis, the same results as would be 
obtained without extended precision. 

The FORTRAN logical operators . eq . . It . . le . . gt . . ge . and their 
analogs in C and Pascal return false if either operand is a NaN. With -f 68 881 
and -ff pa, all except .eq. and . ne . also cause an exception if either 
operand is a NaN. 

Certain features of the IEEE Standard are not readily available from higher-level 
languages. Standard language or library extensions will eventually be developed 
to allow machine-independent access to these features. Sun supports such stan- 
dardization efforts and, in the meantime, provides some provisional assembly 
language methods for accessing the special features. Since these provisional 
methods are nonstandard and probably transitory, limit their use to situations in 
which there is no other way to achieve the intended effect. 
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Rounding Modes and Trap 
Enabling 



Trap Enabling 



Using Special Features From 
FORTRAN 



Rounding modes and trap enabling are available with - f 6 8 8 8 1 and -f f pa. 
These fields are contained in the MC6888 l’s f per register; see the MC68881 
User’s Manual. The following libc routine exchanges a new 32-bit mode con- 
trol word with the existing one: 

as : 

. 

.globl _fpmode_ 

_fpmode_ : 

| On entry, sp0 (4) contains the address of the 
| new 32-bit value to be placed in fper. 

| On exit, dO contains the previous fper value. 

V > 



til : 




cc : 



f 


int 


fpmode_(new) 






int 


*new ; 




s 









pc: 



function fpmode ( new : integer ) : integer ; 

external fortran ; 

s > 



Rounding direction modes are effective with -f 68881 and -f fpa. Rounding 
precision modes are only effective with - f 6 8 8 8 1 . 

Trap enabling with __f pmode_ will cause SIGFPE signals when the 
MC6888 1 detects an exception whose trap is enabled. The SIGFPE code indi- 
cates the specific exception; the MC6888 l’s f piar register contains the 
address of the MC6888 1 instruction which caused the exception (the pc is usu- 
ally advanced beyond that point before the exception is recognized). 



c Set rounding direction mode to round toward zero: 

integer oldmode, fpmode 
oldmode = fpmode (16) 



c Set rounding precision mode to double precision: 

integer oldmode, fpmode 
oldmode = fpmode (12 8) 
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f — \ 

c Enable signalling SIGFPE on overflow, division by zero, 
c or invalid operation: 

integer oldmode, fpmode 
oldmode = fpmode ( 62464) 

V. s 



Exceptions, Condition Code, 
and Quotient Status 



-f 68 8 8 1 and -f f pa provide exceptions, condition code, and quotient status 
in the MC68881’s fpsr register; see the MC6888 1 User’s Manual. The libc 
routine _f pstatus_ exchanges a new 32-bit status word with the existing 
one, somewhat like _f pmode_: 

r — \ 

.globl _fpstatus_ 

_fpstatus__: 

| On entry, sp@(4) contains the address of the 

I new 32-bit value to be placed in fpsr. 

| On exit, dO contains the previous fpsr value. 

\ > 



Only the Accrued Exception byte of fpsr is defined with -ffpa. Here is an 
example of testing the Accrued Exception byte from FORTRAN: 

c Test accrued exception bits to determine if overflow, 
c division by zero, or invalid operation has occurred; 
c also clear status bits: 

integer oldstatus, fpstatus 
oldstatus = fpstatus (0) 

if (and (oldstatus, 128) .ne. 0) print ' invalid operation occurred' 
if (and (oldstatus, 16) .ne. 0) print *,' division by zero occurred' 
if (and (oldstatus, 64) .ne. 0) print *,' overflow occurred' 

\ / 



Signalling NaNs 



Using fpmode and 
f pstatus_ from C 



Signalling NaNs are detected by the MC6888 1 as NaNs containing a cleared 
leading fraction bit When detected, SIGFPE is signalled with the 
F P E_F L TN AN_T RAP code. The FPA always causes recomputation when a NaN 
is encountered; the MC6888 1 distinguishes quiet from signalling during recom- 
putation. -f soft and -f sky treat all NaNs as quiet 

The following program comprehensively demonstrates using the _f pmode_ 
and _£ pstatus_ functions from C to determine existing settings and install 
new ones. 









enum fp rounding direction 
{ 

fp_rd_nearest = 0, 


/* rounding direction type */ 




fp_rd_zero = 1, 
fp_rd__plus =2, 
fp rd minus = 3 
} ; 






V. 




.a 
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/* extended rounding precision type */ 



enum fp_rounding_precision 

{ 

fp_rp_ext ended = 0, 
fp_rp_s ingle = 1 , 
fp_rp_double - 2 , 
fp_rp_3 = 3 

} ; 

enum fp_accrued_type /* accrued exceptions according to fpsr bit number */ 



{ 

fp_inexact = 3, 
fp_divide = 4, 
fp_underflow = 5, 
fp_overflow = 6, 
fp_in valid = 7 
} ; 



enum fp_exception_type /* exceptions according to fpcx/fpsr bit number */ 

{ 



fp__inexl 


= 8 , 


fp_inex2 


= 9 / 


fp_dz 


= 10 


fp_unf 1 


= 11 


fp ovfl 


= 12 


fp__operr 


= 13 


fp snan 


- 14 


fp_bsun 
1 ; 


- 15 



unsigned 

print_81_mode (newmode) /* Exchanges floating-point mode and 

prints out new mode settings. */ 

uns igned newmode ; 

{ 

unsigned oldmode , fpmode__() ; 

oldmode = fpmode__(& newmode) ; 
printf(" New floating-point mode: \n” ); 
printf(" Rounding direction toward ") ; 
switch ( (newmode » 4) & 3) { 

case fp_rd__nearest : 

printf (” nearest ") ; 
break; 

case fp_rd_zero: 

printf (" zero M ) ; 
break; 

case f p_rd_minus : 

printf (" minus infinity ") ; 
break; 

case fp_rd_jplus : 

printf (" plus infinity ") ; 
break; 

} 
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print f (”\n") ; 

printf ( M Extended rounding precision ") ; 
switch ( (newmode » 6) & 3) { 

case f p_rp_extended : 

print f(" extended ") ; 
break; 

case fp_rp__single : 

print f ( M single ”) ; 
break; 

case f p__rp_double : 

print f(" double "in- 
break; 

} 

print f ( M \n”) ; 

printf (" Enabled exceptions: ") ; 
if ( (newmode » (unsigned) fp_inexl) & 1) 
printf ( " inexl " ) ; 

if ( (newmode » (unsigned) fp_inex2) & 1) 
printf (" inex2 ") ; 

if ( (newmode » (unsigned) fp__dz) & 1) 
printf (" dz ") ; 

if ( (newmode » (unsigned) fp_unfl) & 1) 
printf (" unfl ") ; 

if ( (newmode » (unsigned) fp_ovfl) & 1) 
printf (" ovfl ") ; 

if ( (newmode » (unsigned) fp__operr) & 1) 
printf (" operr "); 

if ( (newmode » (unsigned) fp_snan) & 1) 
printf (" snan ” ) ; 

if ((newmode » (unsigned) fp_bsun) & 1) 
printf (" bsun ") ; 
printf (" \n ") ; 
return (oldmode) ; 



unsigned 

print_81_status (newstatus) /* Exchanges floating-point status and 

prints out old status settings. */ 
unsigned newstatus; 

{ 

unsigned oldstatus, fpstatus__() ; 

printf (" Old floating-point status: \n" ); 
oldstatus = fpstatus_(&newstatus) ; 
printf (” Current exceptions: ”); 
if ( (oldstatus » (unsigned) fp_inexl) & 1) 
printf (" inexl ") ; 

if ( (oldstatus » (unsigned) fp_inex2) & 1) 
printf (" inex2 ”); 

if ( (oldstatus » (unsigned) fp_dz) & 1) 
printf (” dz " ); 

if ( (oldstatus » (unsigned) fp__unfl) & 1) 
printf (" unfl M ) ; 
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fp_snan) 

fp_bsun) 



if ((oldstatus » (unsigned) fp_ovfl) & 1) 
printf( M ovfl ") ; 

if ( (oldstatus » (unsigned) fp_operr) & 1) 
printf (” operr ") ; 

if ( (oldstatus » (unsigned) fp_snan) & 1) 
printf( M snan ”) ; 

if ((oldstatus » (unsigned) fp_bsun) & 1) 
printf ( M bsun ,f ) ; 
printf (”\n” ) ; 

printf(" Accrued exceptions: ") ; 
if ((oldstatus » (unsigned) fp_inexact) & 1) 
printf (" inexact ") ; /* inexl + inex2 */ 

if ( (oldstatus » (unsigned) fp_divide) & 1) 
printf (" divide ") ; /* dz */ 
if ( (oldstatus » (unsigned) fp_underf low) & 1) 
printf (” underflow ") ; /* unfl */ 

if ( (oldstatus » (unsigned) fp_overflow) & 1) 
printf (" overflow ") ; /* ovfl */ 

if ( (oldstatus » (unsigned) fp_invalid) & 1) 

printf ( ,f invalid ") ; /* operr + snan + bsun */ 

printf ( f, \n") ; 
return (oldstatus) ; 



main ( ) 

{ 

double 

double 



/* Demonstration program */ 
x; 

zero = 0.0 f one = 1.0, three = 3.0, big = l.Oe+300, 
small = 1.0e-300; 



pr int_8 l_mode ( 0 ) ; 
x = one / three; 



/* Default modes. 
/* Inexact. */ 



print_81_status (0) ; /* Inexact message; clear status. */ 



one / zero; 



/* Divide by zero. */ 



print__81_status (0) ; /* Divide by zero message; clear status. */ 
print_81_mode (0x50) ; /* Round toward zero, round extended to 

* single precision. */ 

x = zero / zero; /* Invalid operation. */ 

x = big * big; /* Overflow, inexact. */ 

x = small * small; /* Underflow, inexact. */ 
print_81_status (0) ; /* Invalid operation, overflow, 

* underflow, inexact message, clear status. */ 



3.3. Special Library Entry 
Points 



Special library entry points are defined to supply certain IEEE operations. The 
names and calling sequences are definitely subject to change, but they are listed 
here to support the IEEE conformance claims below. The first letter indicates the 
type of arithmetic: 

□ V... -f switch 
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□ 


F... 


-f soft 


o 


S... 


-f sky 


□ 


M... 


-f 68881 


□ 


W... 


-ffpa 



while the last letter indicates the precision: 

□ ...s single 

□ ...d double 

Calling conventions are shown in the table below. Not all entry points have two 
arguments. 





First ArgumentIResult 


Second Argument 


Single precision 


dO 


dl 


Double precision 


dO/dl 


a0@ 



The IEEE remainder, a function of two arguments, is available from these assem- 
bly language entry points: 

— 
.globl Vrems, Vremd,Frems, Fremd, Srems, Sremd,Mrems,Mremd, Wrems, Wremd 



Convert floating to integral value in integer format with rounding toward zero, a 

one-argument function, is available from these assembly language entry points: 
. 

. globl Vints, Vintd,Fints, Fintd, Sints, Sintd, Mints, Mintd, Wints, Wintd 

V > 



This function is directly available in FORTRAN (INT), Pascal (trunc), and C 
(int). 



Convert floating to integral value in integer format with rounding to nearest, a 
one-argument function, is available from these assembly language entry points: 




Convert floating to integral value in integer format with rounding according to 
current IEEE rounding mode, a one-argument function, is available from these 
assembly language entry points: 

'i 

• globl Mr int s , Mr intd, Wr int s , Wrintd 

J 



Convert floating to integral value in floating-point format with rounding toward 
zero, a one-argument function, is available from these assembly language entry 
points: 
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.globl Vaints, Vaintd, Faints, Faintd, Saints, Saintd, Maints,Maintd, Waints, Waintd 



This function is directly available in FORTRAN (AINT). 



Convert floating to integral value in floating-point format with rounding to 
nearest, a one-argument function, is available from these assembly language 
entry points: 




Convert floating to integral value in floating-point format with rounding accord- 
ing to current IEEE rounding mode, a one-argument function, is available from 
these assembly language entry points: 

. globl Mar ints , Marintd, Warint s , Warintd 
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Benchmarks 



Several benchmarks are described in this chapter to give you some idea of the 
performance and degree of conformance of the floating-point support provided 
for Sun’s systems. 

The benchmarks described include conformance benchmarks, elementary func- 
tion accuracy benchmarks, and performance benchmarks. 

Note: All specific claims for conformance, accuracy, and speed are based on 
tests and measurements made using a preliminary version of Sun Software 
Release 3.2. Results obtained with the final Release 3.2 or with other 
releases may vary. 

The conformance benchmarks described in this chapter are: 

□ IEEE test vectors. 

□ Paranoia. 

Conformance to the IEEE Standard is primarily tested by running the IEEE test 
vectors, a collection of 15 sets of input data and output numerical results and 
exceptions, distributed by the University of California. A test program converts 
the input data from its distributed symbolic form to the corresponding binary 
representation, obtains the result delivered by Sun floating-point software or 
hardware, and compares the computed result and exceptions with those 
prescribed by the Standard. 

Sun’s test program tests single-precision results computed from single-precision 
operands, and double-precision results computed from double-precision 
operands. The test program obtains access to Sun’s floating-point software or 
hardware primarily through FORTRAN. The following double-precision add test 
function is typical: 

— _ — n 

real* 8 function addd ( x, y) 
real*8 x, y 
addd = x + y 
end 

v J 



Test vectors for the following functions are tested with FORTRAN subprograms 
like addd: 
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□ absolute value 

□ add 

□ compare 

□ divide 

□ multiply 

□ negate 

□ subtract 

□ square root 

Some functions are coded in assembly language as described in Appendix D. 
Test vectors for the following functions are tested with such assembly language 
subroutines: 



□ copysign 

o fraction part 

□ logb 

□ nextafter 

o remainder 

□ round to integer 

□ scalb 



The test vectors as distributed contain 21002 test cases when all rounding modes 
are tested and 5280 test cases when only the default round-to-nearest mode is 
tested. 



Software Test Vector Results 



Sky Test Vector Results 



MC68881 Test Vector Results 



With -f sof t , IEEE test vectors find no numerical errors in the default round- 
ing direction mode. Software floating point does not support other rounding 
direction modes or IEEE exceptions. 

With -f sky, IEEE test vectors find 59 numerical errors in the default rounding 
direction mode, 20 in addition, 20 in subtraction, 8 in multiplication, 2 in divi- 
sion, and 9 in square root. All are errors of one unit in the last place or errors in 
the sign of zero results. The errors arise because of limitations in the microcode 
for the Sky FFP. Sky floating point does not support other rounding direction 
modes or IEEE exceptions. 

With A93N -f 6 8 8 8 1 , and with double rounding precision mode, IEEE test 
vectors find no numerical or exception errors in any rounding direction mode. 

With A93N -f 6 8 8 8 1 , and with extended rounding precision mode, IEEE test 
vectors find 22 numerical errors and 50 exception errors in addition (2, 0), sub- 
traction (2, 0), multiplication (2, 27), division (1 1, 23), and square root (5, 0), 
testing all rounding direction modes. All numerical errors are one unit in the last 
place. The errors arise because extended rounding precision occasionally causes 
different results from double rounding precision as discussed above. 
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A79J MC68881 Test Vector 
Results 



With A79J -f 68 8 8 1, and with double rounding precision mode, IEEE test 
vectors find 27 numerical errors and no exception errors, in multiplication (16), 
division (5), and scalb (6), testing all rounding direction modes. All numerical 
errors relate to erroneous zero results being computed when correctly rounded 
results would be the smallest normalized number. 



Sun FPA Test Vector Results 



Sun FPA Test Vector Results, 

fpamode (3) 



With -f fpa and an A93N MC6888 1 installed, IEEE test vectors find no 
numerical or exception errors in any rounding direction mode. 

With -f fpa and "fast" fpamode enabled, IEEE test vectors find 1616 numerical 
errors and 1032 exception errors in addition, subtraction, multiplication, division, 
comparison, round to integral value, and nextafter. This mode is not intended to 
conform to the IEEE standard. 



4.2. Parano i a Test 
Program 



The Paranoia test program was developed by W. Kahan at the University of 
California to test implementations of floating-point arithmetic by ostensibly 
machine-independent means. Shortcomings in numerics are detected and 
classified as DISASTERS, FAILURES, SEVERE DEFECTS, DEFECTS, 
FLAWS, and WARNINGS. At the end of the program, the number of shortcom- 
ings in each class is displayed, and if the arithmetic appears to conform to the 
IEEE Standard, the program so states. Sun tests Fortran, C, and Pascal versions 
of Paranoia that evolved from Kahan’s original 1983 version in Basic. 



In the following tables, DEFECTS, FLAWS, and WARNINGS encountered are 
enumerated Di, Fj, and Wk, and described afterward. 754 indicates the arith- 
metic appeared to conform to the IEEE Standard. 















Single-Precision Paranoia 


Results 




Compiler Options 


Fortran 


C 


Pascal 




-0 -ffpa 


754 W1 


754 W1 


754 W1 




-0 -f 688 81 




754 W1 


754 W1 




-0 -f68881 -fstore 


754 W1 








“0 -fsky 


Dl D2 D3 FI F2 W1 


D4 FI F2 W1 


Dl FI F2 W1 




-0 -fsoft 


754 W1 


754 W1 


754 W1 




V 








J 



r 




















Double-Precision Paranoia 


Results 




Compiler Options 


Fortran 


c 


Pascal 




-0 


-f fpa 




754 W1 


754 Wl 


754 Wl 




-0 


-f 68881 






Fl Wl 


Fl Wl 




-0 


-f 68881 


-fstore 


Fl Wl 








-0 


-fsky 




FI Wl 


Fl Wl 


Fl Wl 




-0 


-fsoft 




754 Wl 


754 Wl 


754 Wl 






Key to 


Shortcomings : 








Dl 


: DEFECT: 


sqrt ( 0 


2401000e+04) - 0. 


4900000e+02 


0 . 3814697e-05 






instead 


of correct value 0 . 








D2 

< 


: DEFECT: 


computed 


( 0.20000000e+01) 


**(11) = 0 


. 20480002e+04 
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compares unequal to correct 0 . 20480000e+04 ; 
they differ by 0 . 24414063e-03 . 

Error like this may invalidate financial 
calculations involving interest rates. 

D3: DEFECT: computed ( 0 . 20000000e+01) ** ( 120 ) - 0 . 132 92334e+37 

compares unequal to correct 0 . 13292280e+37 ; 
they differ by 0 . 53875151e+31 . 

D4: DEFECT: Comparison alleges that what prints as Z = 1 . 40129846432481710e“45 

is too far from sqrt (Z) ** 2 = 0 . 00000000000000000e+00 . 

FI: FLAW: lack(s) of guard digits or failure (s) to correctly round or chop 
F2 : FLAW: underflow can stick at an allegedly positive value zO 
that prints out as 0 . 14012985e-44 
Wl: WARNING: computed ( 0 . 00000000e+00) ** ( 0) - NaN 

compares unequal to correct 0 . 10000000e+01 

v / 



-f fpa Comments 



That peculiar results for x* *y count only as WARNINGS when x=0 and 
y<=0 reflects a division of opinion among analysts. The original BASIC version 
of Paranoia reflects Kahan’s view that x**0 should be 1, even if x be zero, 
infinity, or NaN. Subsequent translators felt strongly differently, and so changed 
the code accordingly to note these cases as WARNINGS without otherwise 
counting them in determining the quality of implementations. 

In "fast” mode, selected by 



call fpamode(3) 

\ > 



Paranoia discovers an inconsistency of underflow treatment: underflow itself 
leads to recomputation of a subnormal result, but in subsequent use that result is 
sometimes treated as zero. Thus in single precision, Paranoia output 
includes the following: 

SERIOUS DEFECT: X = 1.61630e-38 is Unequal to Z = 1.17549e-38, 

yet X-Z yields 4.40810e-39. 

V, / 



-f 68 881 Comments 



"Fast" mode is not intended to conform to the IEEE Standard. 

The FORTRAN version of Paranoia must be compiled with -f store to 
force rounding to storage precision on assignment. If -f store is omitted. 
Paranoia encounters numerous FAILURES and SERIOUS DEFECTS. 

Double precision does not conform to the IEEE standard even with -f store 
because of multiple roundings. However, by selecting double rounding precision 
mode with 



( — 




a 




call fpmode(128) 




L_ 




j 



the result becomes 754 FI. 
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-f sky Comments 



Shortcomings, particularly in single precision, are due to limited microcode 
space on the Sky FFP. C, FORTRAN, and Pascal defects differ because sqrt 
and pow are implemented in double precision in C, while pow is implemented 
in double precision in Pascal. 



Summaries 



Each program prints an overall summary of the tested arithmetic. The sum- 
maries corresponding to die tables above are as follows: 





754 Wl: 

No failures, defects nor flaws have been discovered. 

Rounding appears to conform to IEEE standard p754. 

The arithmetic diagnosed appears to be excellent ! 

FI Wl: 

The arithmetic diagnosed seems satisfactory though flawed. 

D1 D2 D3 FI F2 Wl : 

D4 FI F2 Wl: 

D1 FI F2 Wl: 

The arithmetic diagnosed may be acceptable despite inconvenient defects. 

J 



4.3. Elementary Function 
Accuracy Benchmarks 



The accuracy of elementary transcendental functions is tested using the 
Elef unt FORTRAN test suite of Cody and Waite, and a preliminary C test suite 
of Alex Liu, of the University of California. More information on these pro- 
grams may be obtained from sources listed in the Preface to this manual. 

All programs are compiled with -0; MC6888 l’s are used only in the default 
extended rounding precision mode. For these tests, A79J and A93N masks pro- 
duce the same results. 



Elef unt Test Programs 



In Cody and Waite’s program, each function is tested by checking several 
independent identities at many random points; reported errors measure the degree 
to which the computed function values fail to satisfy certain identities, rather 
than the errors in the computed function values themselves. In addition to testing 
identities at random points, as shown below, Cody and Waite’s programs also 
test error responses at isolated points. 

Cody and Waite advise that it may be necessary to modify their test programs 
when extended-precision expression evaluation is performed, as in the case of 
Sun’s -f 68881 option. Sun’s version of Elefunt accordingly contains 
changes to the codes that purify test arguments; these changes are listed below 
under the commented-out source lines they replace: 
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c x - (x + eight) - eight alog0870 

x = x + eight 
x = x - eight 

c x = (x + eight) - eight alog0940 

x = x + eight 
x = x - eight 

c x - (x + w) — w powr0940 

X = x + w 
X = x - w 

c y = (x + y) - x sin00760 

y = y + x 
y = y - x 



Below are worst-case differences detected in identities tested, expressed as the 
number of bits of precision lost. 



r 

Greatest number 


of bits lost 


- single- 


-precision 


elefunt 




compiler option-> 


-fsoft 


-f sky 


-f 68881 


-f 68881 
-fstore 


-ffpa 


function : 


sqrt 


0 


1.0 


0 


0 


0 


exp 


1.0 


2.0 


0 


1.0 


1.0 


log 


1.0 


2.6 


0 


1.0 


1.0 


loglO 


2.4 


3.0 


3.0 


2.1 


2.1 


x**y 


1.0 


7.3 


3.2 


1.0 


1.0 


sinh/cosh 


2.0 


2.0 


0 


1.0 


1.6 


tanh 


2.1 


2.0 


0 


1.0 


1.5 


sin/cos 


1.5 


3.1 


0 


0.6 


1.5 


tan 


1.8 


2.6 


0 


1.3 


1.8 


asin/acos 


1.9 


2.0 


0 


1.0 


1.0 


atan 

l 


1.0 


1.9 


0 


1.0 


1.0 

/ 
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/ 

Greatest number 


of bits lost 


- double- 


-precision 


elefunt 




- N 


compiler option-> 


-f soft 


-fsky 


-f 68881 


-f 68881 
-f store 


-ffpa 




function : 


sqrt 


0 


0 


0 


0 


0 




exp 


1.0 


1.0 


0 


1.0 


1.0 




log 


1.0 


1.0 


0 


1.0 


1.0 




loglO 


2.4 


2.6 


3.1 


2.1 


2.1 




x**y 


2.3 


2.4 


6.2 


2.2 


2.2 




sinh/cosh 


2.1 


2.1 


0 


1.0 


1.6 




tanh 


2.1 


2.1 


0 


1.0 


1.5 




sin/cos 


1.5 


1.5 


0 


0.7 


1.5 




tan 


2.0 


2.1 


0 


1.2 


2.0 




asin/acos 


1.2 


1.2 


0 


1.0 


1.0 




atan 

V 


1.0 


1.4 


0 


1.0 


1.7 


J 



Test Programs of Liu 



In Liu’s test program, 64 random points are tested in 64 regions for each func- 
tion; reported errors are the errors in the computed function values themselves, 
tested against more accurate values generated internally by Liu’s program. Tests 
of e**x-l and log(l+x) used assembly language entry points 
[FSMW]{exp 1 ,log 1 } [sd]. 



Below are worst-case errors in ulps, units in the last place, detected for each 
function tested. A difference of one ulp corresponds to a difference in the least 
significant bit. 



— 
Greatest errors 


in ulps - 


single-precision 


Liu tests 




1 \ 


compiler option-> 


-fsoft 


-fsky 


f 68881 


■ffpa 




function : 


sin 


0.88 


5.79 


0.54 


0.88 




cos 


0.90 


5.95 


0.53 


0.90 




atan 


0.86 


2.08 


0.52 


0.86 




log 


0.92 


3.48 


0.53 


0.92 




log (1+x) 


0.88 


0.90 


0.53 


1.06 




e**x - 1 
v 


0.94 


0.94 


0.53 


0.95 


J 
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f 

Greatest errors 


in ulps - 


double-precis ion 


Liu tests 






compiler option-> 


-fsoft 


-fsky 


f 68881 


-f fpa 




function : 


sin 


1.01 


0.99 


0.57 


1.01 




cos 


1.01 


0.90 


0.55 


1.01 




atan 


1.06 


0.98 


0.56 


1.02 




log 


0.92 


0.96 


0.59 


0.92 




log (1+x) 


0.83 


0.83 


0.59 


0.83 




e**x - 1 

V 


1.13 


1.00 


0.57 


1.18 


j 



Monotonicity 



Operations specified in the IEEE Standard ate required to be monotonic; elemen- 
tary transcendental functions are often expected to be monotonic in their primary 
domains. Although expectations about monotonicity are seldom explicit in pro- 
grams, monotonicity failures can befuddle debugging. 



Liu’s program tests monotonicity at selected points; failures were detected only 
with -fsky and only in single precision as follows: 



f 

log: 


4 monotonicity 


failures 


including 


3 . 2768000000000000e+04 




atan: 

V 


449 monotonicity 


failures 


including 


5 . 222 527 8 12 50 00 00 0e+05 


j 



Although exhaustive testing of monotonicity is infeasible for double precision, it 
is in principle possible for single precision. At present, the following single pre- 
cision IEEE Standard and Fortran functions have been tested over the indicated 
primary domains by evaluation at each single-precision representable number in 
the domain. If monotonic, execution times are shown in hours; "no" means not 
monotonic: 

— 



IEEE function 


and domain 






-f fpa 


-f 68881 


-fsky 


-fsoft 














A93N 






x*x 


[ 


0. 


Inf 


] 




8 


30 


25 


abs 


[ 


0. 


Inf 


] 


3 


5 


17 


7 


aint 


[ -Inf 


r 


Inf 


3 


15 


14 




20 


arint 


[ -Inf 


r 


Inf 


3 


20 


34 






ceil 


[ -Inf 


t 


Inf 


3 


22 


46 






floor 


[ -Inf 


r 


Inf 


3 


21 


46 






dble 


[ -Inf 


r 


Inf 


3 


14 


10 


72 


37 


int 


[ -0 . 2147 4836e+10 , 


2 . 14748352e+9 


3 




18 






rint 


[ -0.21474836e+10, 


2 . 14748352e+9 


3 




24 






sqrt 


[ 


o. 


Inf 


3 


10 


9 


71 


49 
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FORTRAN function and domain 



-ffpa -f 68881 -f sky -fsoft 
A93N 



acos 


[ -1.0000000 


f 


1.0000000 


i 


29 


no 


no 


anint 


[ -Inf 


r 


Inf 


i 


28 


43 




asin 


[ -1.0000000 


r 


1.0000000 


i 


21 


19 


no 


atan 


[ -Inf 


f 


Inf 


i 


16 


26 


no 


cos 


[ -3.1415925 


r 


0. 


i 


4 


9 


72 


cosh 


[ o. 


r 


Inf 


i 




22 




exp 


[ -Inf 


t 


Inf 


i 




41 


no 


exp-1 


[ -Inf 


r 


Inf 


i 




62 




2**x 


[ -Inf 


r 


Inf 


i 




36 




10**x 


[ -Inf 


r 


Inf 


i 




44 




log 


[ o. 


f 


Inf 


i 


19 


31 


no 


logl+x 


[ -1. 


r 


Inf 


i 


22 


46 


no 


log2 


[ o. 


r 


Inf 


i 




38 




loglO 


t o. 


f 


Inf 


i 




27 




nint 


[ -0 .2147 483 6e+10 


f 


2 . 14748352e+9] 


27 


27 




sin 


[ -1.5707963 


f 


1.5707963 


i 


4 


11 


no 


sinh 


[ -Inf 


r 


Inf 


i 




44 




tan 


[ -1.5707963 


r 


1.5707963 


i 


16 


14 


90 


tanh 


[ -Inf 


r 


Inf 


i 
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4.4. Performance 
Benchmarks 

Linpack 



The -f 68881 acos monotonicity failure is in extended precision and is not 
detected if -f store is in effect. 

All tests were compiled with -0 and without -f store, and were made on 
Sun-3’s except -f sky, made on a Sun- 2. Timings are the entire execution time 
for the test program and consequently give only a rough comparison of function 
speed. As die execution times imply, these tests constitute an ongoing project. 

The performance benchmarks described in this section, Linpack and SPICE, 
are commonly used to compare the relative performance of different 
manufacturer’s computers. 

The most realistic of the frequently cited benchmarks of peak floating-point 
speed is the Linpack benchmark, a FORTRAN program that determines the 
time required to solve a 100x100 system of linear equations using the Linpack 
linear algebra subroutine library. 

Sun’s version of the Linpack benchmark program reports performance in 
thousands of floating-point operations per second (KFLOPS). Performance 
reported below is based upon measuring total user and system CPU time used. 

All results are based on compiling FORTRAN source code; -fsoft, -fsky, 
and -f switch were compiled on Sun-2’s; -f68881 and -ffpa were com- 
piled on Sun-3 ’s. 

All floating-point hardware runs at 16.7 MHz except as noted. 
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Sun-2 results (10 MHz MC68010 CPU): 

Single Double 

Compiler Precision Precision 

Options KFLOPS KFLOPS Comments 

-O -fsoft 13 6 

-0 -fswitch 32 18 using Sky FFP 

-O -fsky 49 27 



Sun-3 results (15 MHz MC68020 CPU): 



r 




Single 


Double 




FPU 


Compiler 


Precision 


Precision 




MHz 


Options 


KFLOPS 


KFLOPS 


Comments 




-0 -fsoft 


28 


12 




12.5 


-0 -fswitch 


47 


34 


using MC68881 


12.5 


-0 -f 68 881 


78 


73 






-0 -fswitch 


55 


38 


using MC68881 


16.7 

V 


-0 -f 68881 


93 


87 


j 



Sun-3 results (16.7 MHz MC68020 CPU): 









Single 


Double 




FPU 


Compiler 


Precision 


Precision 




MHz 


Options 


KFLOPS 


KFLOPS 


Comments 




-0 


-fsoft 


34 


15 




12.5 


-0 


-fswitch 


50 


37 


using MC68881 




-0 


-fswitch 


65 


45 


using MC68881 


12.5 


-0 


-f 68881 


85 


80 




16.7 


-0 


-f 68881 


108 


100 






-0 


-fswitch 


185 


105 


using FPA 




-0 


-ffpa 


500 


315 




16.7 


-0 


-f 68881 


87 


82 


double rounding 


16.7 


-f 68881 


95 


87 




16.7 


-0 


-f 68881 










-fstore 


99 


90 




16.7 


-0 


-f 68881 


109 


101 


rolled 




-f fpa 


280 


160 






-P 


-ffpa 


420 


285 






-0 


-ffpa 


615 


405 


rolled 



The comment "double rounding" means that the double rounding precision mode 
of the MC6888 1 was selected. The comment "rolled" means that the inner loops 
of the Basic Linear Algebra Subroutines were rewritten slightly. The distributed 
form of the Linpack benchmark program has the key inner loop written in the 
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"unrolled" form 



do 1 i = 1, n, 4 
x (i ) = x (i ) + c 
x(i+l) = x(i+l) + c 
x(i+2) = x(i+2) + c 
1 x(i+3) = x(i+3) + c 



* y (i ) 

* y(i+l) 

* y(i+2) 

* y(i+3) 



because that form was faster on certain mainframes common in the mid-1970’s. 
However, the unrolling defeats many current vectorizing compilers, so super- 
computer manufacturers usually measure the rolled speed by rewriting the loop 
in the form: 



f 


do 1 


i = 1, n 


\ 


1 


x (i 


) = x (i ) + c * y (i ) 










> 



Dongarra’s Linpack benchmark performance compilation, listed in the Pre- 
face, annotates such results as "Rolled BLAS." Further complicating the issue is 
that compilers on some systems do not generate optimum code for the inner loop 
whether rolled or unrolled, so hand-coded assembly language is faster yet; these 
results are annotated as "Coded BLAS.” For the rolled form of the Linpack 
benchmark, the Sun-3 code generated with -O and either -f f pa or -f68881 
is truly optimal and cannot be improved by hand coding in assembly language. 

The SP ICE program simulates integrated circuits and is an example of a com- 
plete application that is floating-point intensive but does not continually attain 
peak performance rates as does the Linpack benchmark. It is therefore more 
representative of relative performance on realistic computations that do not 
closely fit the model of a Linear Algebra computation. 

For benchmarking purposes Sun uses a FORTRAN implementation of SPICE 
version 2G.6, and a C implementation of SP ICE version 3A.7, both with an 
input data file called MOSAMP2, listed in Appendix E. The reported time is the 
total elapsed real time in seconds for the complete program, invoked with a com- 
mand like 

' > 
time spice. out < mosamp2 . input >! mo s amp2 . output 
>. 



In each case the program was run at least twice and the fastest result reported. 
Note that the global optimizer iropt is available in Sun’s f77,butnotin 
Sun’s cc. -f switch and -f soft versions were compiled on Sun-2’s and 
run on Sun-2s and Sun-3s. -f sky versions were compiled and run on Sun-2s. 
- f 6 8 8 8 1 and - f f pa versions were compiled and run on Sun-3s. All 
floating-point hardware runs at 16.7 MHz except as noted. 
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Results are: 








r 






\ 


FPU Compiler Real 

MHz Options 


Time to 
(Seconds) 


Execute 




Sun-2 results (10 MHz MC68010 CPU) : 


2G . 6 
Time 


3A.7 

Time 




-0 -fsoft 


1372 


1979 




— 0 -fswitch+Sky 


483 


674 




-0 -fsky 


441 


605 




Sun-3 results (15 MHz MC68020 CPU) : 


-0 -fsoft 


595 


858 




12.5 -0 — f switch+A7 9 J 


184 


262 




12.5 -0 -f 68881+A7 9 J 


102 


128 




16.7 -0 — f switch+A93N 


166 


237 




16.7 -0 -f 68881+A93N 


88 


109 




Sun-3 results (16.7 MHz MC68020 CPU): 


-0 -fsoft 


482 


696 




12.5 -0 -fswitch+A79J 


165 


233 




12.5 -0 -f 68881+A79J 


92 


117 




16.7 -0 -f switch+A93N 


139 


198 




16.7 -0 -f 68881+A93N 


78 


96 




-0 -fswitch+fpa 


76 


105 




-0 -ffpa 

V 


45 


60 


j 



4.5. Benchmarking 
Hazards 



MC68020 Cache 



It is not easy to construct meaningful benchmark programs that are short and 
easy to understand and provide a sound basis for projecting results for realistic 
applications. The following examples illustrate how seemingly minor variations 
in hardware or in coding techniques can have surprising results. 

The MC68020 has a 256-byte direct-mapped instruction cache. "Direct-mapped" 
means that memory address X is mapped to cache address X mod 256, so 
memory locations X and X+2 5 6 * n cannot be in the cache at the same time. 
This can be a problem for short loops or subroutines that call other short subrou- 
tines. If the difference in memory addresses of the caller and callee happens to be 
close to a multiple of 256, almost every instruction fetch misses the cache and 
performance is degraded overall. Note that this is a problem particularly of 
direct-mapped caches, and is more likely to be encountered when such caches are 
small, and is most noticeable when the caller and callee are very short, so that the 
instruction fetches from main memory add up to a significant fraction of the 
overall time required. 

Programs coded in higher-level languages or even in separate assembly-language 
modules are usually written without thought to where they might be placed in 
memory, since that is often an accident of code generation or linking order which 
may not be visible to or controllable by the programmer. 
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The Whetstone benchmark was originally an Algol 60 program that dupli- 
cated the typical instruction stream processed by the Whetstone Algol interpreter 
during the 1960’s. Performance is reported as thousands of Whetstone inter- 
preter instructions per second. 

Various FORTRAN translations of this program have been used to benchmark 
floating-point performance, but interpretation of the results has become increas- 
ingly difficult as new computer architectures and optimization techniques 
develop. The Linpack and SPICE benchmarks discussed above give better 
information about relative performance of computers and compilers. However, 
the Whetstone program is ideal for illustrating MC68020 cache effects. 

The Whet stone program consists of several timing loops, but the most impor- 
tant, that consumes about half the overall time, is a loop that calls a subroutine 



P3. The loop is 


r 


N 


DO 90 1=1, N8 
CALL P3 (X, Y, Z) 

90 CONTINUE 




> 


and the subroutine is 




\ 


SUBROUTINE P3(X,Y,Z) 

COMMON T,T2 
XI = X 

Y1 = Y 

XI = T * (XI + Yl) 

Y1 = T * (XI + Yl) 

Z = (XI + Yl) / T2 

RETURN 

END 

^ ... ... 


-J 



This code is intended to model an important class of numerical codes that solve 
nonlinear problems: optimization and finding zeros of functions. These codes 
repeatedly call user-supplied subroutines to evaluate the nonlinear functions in 
questions. Usually, the calling routine is rather more complicated than a DO 
loop and the parameters vary; in those respects P3 is a deficient abstraction. 

When compiled in single precision with -O and -f f pa, the results obtained 
for this program vary markedly depending on the cache contention between P3 
and its calling loop. Single precision and -f f pa were chosen for illustration 
since cache contention is most noticeable when the timings of the caller and cal- 
lee are as similar as possible. Worst-case performance is 2100 thousand Whet- 
stone instructions per second (KWIPS), best case is 2400. Which result is 
obtained depends on the relative compiled addresses of the do loop and P3, 
which vary among compiler releases. 

The performance with the MC68020’s cache disabled is about 1750 KWIPS, so 
even the worst cache case is much better. A more complicated cache architecture 
on the MC68020 might show less variation but might be slower overall. 
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Assembly Language Inline 
Expansion 



One way to avoid cache contention between callers and callees is to expand pro- 
cedures inline in the calling code. There is a facility in Sun’s compilers to facili- 
tate such expansion at the assembly-language level; see the appendix 
"Assembly-Level Inline Expansion". This facility is intended to allow short, 
simple assembly-language library routines to be expanded inline to bypass the 
procedure-call overhead of saving and restoring registers and pushing arguments. 
To demonstrate the facility, we apply it to P 3 by first compiling P 3 separately 

with -0 -ffpa -S to obtain this assembly language: 



.data 

. comm BLNK , 4 0 

.text 

.globl _j>3_ 

link a6,#-16 
fpmoved fpa4, a6@ (-16) 
fpmoved fpa5,a6@(-8) 
movl a 60 ( 8 ) , aO 
fpmoves aO0,fpa4 
movl a6@ (12) , aO 
fpmoves aO0,fpa5 

f pams f pa 5 , f pa 4 , BLNK , f pa 4 

fpams fpa5,fpa4, BLNK , fpa5 

fpadd3s fpa5, fpa4, fpa2 

fpdivs BLNK +8 , f pa2 

movl a60 (16) , aO 
fpmoves fpa2,a06 
fpmoved a60(-8),fpa5 
fpmoved a60 (-16) , fpa4 
unlk a 6 
rts 

< J 



which was converted manually according to the inline expansion rules to 



.data 




. comm 


BLNK , 40 


.text 




. inline 


_P 3_, 12 


movl 


sp0 + , aO 


fpmoves 


a06, fpaO 


movl 


sp0+, aO 


fpmoves 


aO0 , fpal 


fpams 


fpal, fpaO, BLNK 


fpams 


f pa 1 , f pa 0 , BLNK 


fpadd3s 


fpal, fpaO, fpa2 


fpdivs 


BLNK +8, fpa2 


movl 


sp0+, aO 


fpmoves 


fpa2, aO0 


.end 





fpaO 

fpal 
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When the FORTRAN Whetstone source, with P3 removed, was compiled 
with a . il inline expansion file containing the above, the resulting performance 
was 3300 KWIPS. The inner loop was: 



r 








L77069: 








movl 


#VAR_SEGl+4,aO 






movl 


d6, BLNK +36 






movl 


d7, BLNK +32 






movl 


d3, BLNK +28 






fpmoves 


aO0, fpaO 






movl 


#VAR SEG1+8, a0 






fpmoves 


aO0, fpal 






fpams 


fpal, fpaO, BLNK 


, fpaO 




fpams 


f pa 1 , f pa 0 , BLNK 


_, fpal 




fpadd3s 


fpal, fpaO, fpa2 






fpdivs 


BLNK +8, fpa2 






movl 


#VAR_SEGl+12,aO 






fpmoves 


fpa2,a0@ 






movl 


BLNK +32, d7 






movl 


BLNK +36, d6 






movl 


BLNK +28 , d3 






dbra 


d5, L77069 






s 






J 



Source Level Inline Expansion 



Of course, someone optimizing by hand would normally expand FORTRAN 
source inline into the FORTRAN source of the caller! If the following source 
code is compiled: 

. 

DO 90 1=1, N8 
XI = X 

Y1 = Y 

XI = T * (XI + Yl) 

Y1 = T * (XI + Yl) 

Z = (XI + Yl) / T2 
90 CONTINUE 



then resulting performance is 4100 KWIPS. Examination of the generated code 
shows the same number of floating-point instructions executed: 







A 


L77069: 


fpams05 


fpa5, fpa7, BLNK , fpa6 




fpams05 


fpa5,fpa6, BLNK , fpa4 




fpadd3s05 


fpa4, fpa6, fpa8 




fpdivs 05 


BLNK +8, fpa8 




dbra 


d5, L77069 




v 




> 
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Global Interprocedural 
Analysis 



Performance, Source Coding, 
and Optimization 



Interestingly enough, executing the P3 loop only once would be logically 
sufficient. Inspection shows that the inputs X, Y, T, and T2 are the same on 
each invocation, so the call could be moved outside the loop entirely. An optim- 
izing compiler that did interprocedural global analysis might detect this. To 
measure the effect such analysis would have if performed by a hypothetical com- 
piler, we change the source code to: 



f 








CALL P3 (X, Y, Z) 






DO 90 1=1, N8 




90 


CONTINUE 








J 



Resulting performance is 4800 KWIPS; the floating point instructions are only 
executed once. 

There is no compiler that can optimize code as well as a person who understands 
what is required to be done and how to do it! The Whetstone benchmark exam- 
ple illustrates this point extremely well, simultaneously rendering suspect its own 
value as a general-purpose floating-point benchmark. 

A benchmark program by Myron Ginsberg was adapted to further illustrate how 
source coding techniques and compiler optimization techniques interact to affect 
the throughput on floating-point intensive calculations. This adaptation measures 
the time to evaluate a specific fifth-degree polynomial A at 500 points X ( I ) 
using several different methods. The "explicit" method is written out in extenso 

specifically for fifth-degree polynomials: 
— 

DO 10 1=1,500 

10 Y (I) = A ( 1 ) + X(I) * (A(2) + X(I) * (A (3) + X(I) * 
(A<4) + X(I) * (A(5) + X (I) * A(6) ) ) ) ) 
s i 



The "external" method relies on a conventional external subroutine for evaluating 
polynomials of arbitrary degree: 

> 

DO 10 1=1,500 
Y (I) = PEVAL ( A, M, X(I) ) 

10 CONTINUE 

REAL FUNCTION PEVAL ( A, M, X ) 

REAL A(*) ,X 
INTEGER M,K 
REAL P 
P = A (M+l ) 

DO 15 K = M, 1,-1 
P = P * X + A(K) 

15 CONTINUE 
PEVAL = P 
END 

» j 



The "inline" method effectively expands the "external" code inline at the source 
level: 
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( 












P = A (NCOEF) 








DO 10 1=1,500 








DO 20 J=M, 1,-1 






20 


P = A (NCOEF— J) + X(I) * P 






10 


Y(I) = P 




s 






J 



The "vectored" method is intended to be optimized into appropriate vector opera- 
tions on supercomputers: 



DO 10 L=l,500 

10 BR (L, NCOEF) = A(NCOEF) 

DO 20 J=1,M 
DO 20 1=1,500 

20 BR ( I , NCOEF- J) = A(NCOEF-J) + X(I) * BR(I,NCOEF-J+l) 

DO 30 LL=1, 500 
30 Y (LL) = BR (LL, 1 ) 



The results are reported in KFLOPS. Each polynomial evaluation requires ten 
floating-point operations. The number of floating-point operations executed is 
the same in every case: 





Single-Precision KFLOPS 


\ 


Options 


explicit 


inline 


external 


vectored 


-O -f sky 


62 


49 


43 


23 


-O -f 68881 


210 


130 


125 


73 


-O -ffpa 


1500 


710 


295-335 


175 

j 



r 


Double-Precision KFLOPS 


a 


Options 


explicit 


inline 


external 


vectored 


-O -fsky 


33 


26 


24 


13 


-O -f 68881 


200 


120 


115 


67 


-O -ffpa 
- 


1060 


480 


265-290 


130 

J 



Note that the MC68020’s internal cache causes the variable external case perfor- 
mance with -ffpa. The variation is insignificant with -f 68881. 

As the results indicate, the Sun FPA attains almost one third of the maximum 
theoretical throughput of the 1164/1165 chips - 4800 KFLOPS single, 3300 
KFLOPS double. But this remarkably high efficiency is seldom achievable on 
realistic problems because performance is primarily limited by bus and memory 
bandwidth. The explicitly coded polynomial evaluation creates expressions so 
simple that Sun’s FORTRAN compiler can allocate all the constants to FPA regis- 
ters to minimize bus traffic. The vectorization-oriented evaluation, in contrast, is 
more difficult to optimize effectively for the FPA, since its complicated 
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expressions involve double subscripts and heavy bus traffic. These differences 
are most pronounced for high-performance floating-point hardware. 
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adb Changes 



A.l. Changes in Release 3.1 Release of the Floating-Point Accelerator (FPA) required some changes to adb, 

in order to support assembly language debugging of programs that use the FPA. 
Here are the changes made to adb in Release 3.1: 

1. The new debugger variables A through Z are reserved for special use by 
adb. They should not be used in adb scripts. 

2. The FPA register names f paO through f pa31 are recognized and can be 
used or modified in debugger commands. (This extension only applies to 
machines with FPA’s.) 

3. The debugger variable F governs FPA disassembly. A value of 0 indicates 
that all FPA instructions are to be treated as move instructions. A nonzero 
value indicates that FPA instruction sequences are to be disassembled and 
single-stepped using FPA assembly language mnemonics. The default value 
is 1 on machines with FPA’s; on other machines, the default value is 0. 

4. The debugger variable B is used to designate an FPA base register. If FPA 
disassembly is disabled (the F flag = 0) its value is ignored. Otherwise, its 
value is interpreted as follows: 

0 through 7: 

Based-mode FPA instructions that use the corresponding address regis- 
ter in [ aO . . a7 ] to address the FPA are also disassembled using FPA 
assembler mnemonics. Note that this is independent of the actual run- 
time value of the register. 

otherwise: 

All based-mode FPA instructions are disassembled and single-stepped 
as move instructions. 

The default value of the FPA base register number is —1, which designates 
no FPA base register. 

5. The command $x has been added to display the values of FPA registers 
fpaO through f pa 15, along with FPA control registers and the current con- 
tents of the FPA instruction pipeline. All registers are displayed in the for- 
mat: 

''l 

<low word> chigh word> <double precision> <single precision> 

V / 
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A.2. Examples of FPA 
Disassembly 



This verbose display is used because FPA registers are typeless; in particular, 
they may contain either single- or double-precision floating-point values. If a 
single-precision value is stored, it is always stored in the high-order word. 
Machines without an FPA display the message “no FPA”. 

1. The command $X is similar to $x, but displays the FPA registers fpal 6 
through fpa31 instead of fpaO through fpal5. This is done as a separate 
command because adb cannot display the contents of all FPA registers in a 
single standard-size window. 

2. The command $R displays the contents of the data and control registers of 
the standard MC6888 1 floating-point coprocessor. 



As an example, consider the following assembly source fragment: 




On machines without an FPA, the default mode is to disassemble all FPA 
instructions as moves. For the example program, the following output is pro- 
duced (except the parenthesized comments added here for explanation): 

% as foo.s -o foo.o 

% adb foo.o 

<F=d 

0 (default value of ’F’ on a machine without FPA) 

foo?ia 

foo: movl dO, 0xe0000380 (normal disassembly) 

s, > 



FPA disassembly can be enabled by setting the debugger variable F to 1. For 
example: 



r 

% adb foo.o 

1>F 

<F=d 

1 (new value of ’F’) 

foo?ia 




A 


foo: fpadds dO,fpaO 


(FPA disassembly) 








-j 



On machines with an FPA, FPA disassembly is on by default, so the above out- 
put is produced without having to set the value of F. 

FPA instructions may address the FPA using a base register in [ aO . . a7 ] . In 
practice, only [ aO . . a 5 ] are used by the compilers. 

adb does not know which register (if any) is being used to address the FPA in a 
given sequence of machine code. However, another debugger variable (B) may 
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be set by the user to designate a register as an FPA base register. By default, this 
variable has the value — 1, which means that no register should be assumed to 
point to the FPA, so only instructions that access the FPA using absolute address- 
ing are recognized as FPA instructions. 



For the example program, a machine with an FPA produces the following output: 



r 

% adb foo.o 
<F=d 






\ 


1 

<B=d 

-1 

foo, 3?ia 


(default value of ’F* on a machine with FPA) 
(default value of ’B’) 




foo: 


fpadds dO/fpaO 


(FPA disassembly) 




0x6 : 


movl dO, a0@ (0x380) 


(normal disassembly) 




Oxa: 

Oxe: 

^ _ 


movl dO, a5@ (0x380) 


(normal disassembly) 


J 



Note that the second and third instructions are still disassembled as moves, since 
adb cannot assume that they access the FPA. Continuing this example, if the 
FPA base register number is set to 5, the following output is produced: 



/ 

% adb foo.o 

5>B 

<B=d 

c 






A 


foo, 3?ia 


foo : 


fpadds d0,fpa0 


(FPA disassembly) 




0x6 : 


movl d0,a0@ (0x380) 


(normal disassembly) 




Oxa : 
Oxe : 


fpadds@5 dO/fpaO 


(FPA disassembly) 





Note that the second instruction is still disassembled as a move, since a5, the 
register designated as the FPA base, is not used in it. 

A.3. Examples of FPA FPA data registers can be displayed using a syntax similar to that used for the 

Register Use MC6888 1 coprocessor registers. Note that unlike the MC6888 1 registers, FPA 

registers may contain either single-precision (32-bit) or double-precision (64-bit) 
values; MC68881 registers always contain extended-precision (96-bit) values. 



For example, if f paO contains the value 2.718282, we may display it as follows: 



r 








A 


<fpa0=f 


fpaO 


0x402df 855 


+2.718282e+00 




c. . ... 








V 



Note that the value is displayed in hexadecimal as well as in floating-point nota- 
tion. Unfortunately, an FPA register can only be set to a hexadecimal value. To 
set f paO to 1.0, for example, you must know that this is represented as 
0x3f800000 in IEEE single-precision format: 
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dbx and dbxtool Changes 



B.l. Changes in Release 3.1 Release of the Floating-Point Accelerator (FPA) required some changes to dbx 

and dbxtool, in order to support debugging of programs using the FPA. Here 
are the changes made to dbx in Release 3.1: 

1. There is a new f paasm debugger variable to control disassembly of FPA 
instructions. This variable may be set or displayed using the dbxenv com- 
mand, for which the syntax is: 

f \ 

dbxenv fpaasm [on | off] 

v J 



All FPA instructions are disassembled as moves if the value of fpaasm is of f . 
FPA instructions are disassembled with FPA assembly language mnemonics if 
the value is on. Defaults: on a machine with an FPA, fpaasm is initially set to 
on; on machines without an FPA, it is initially set to of f . 

1. The f pabase debugger variable has been added. It designates an 

MC68020 address register for FPA instructions that use base+short displace- 
ment addressing to address the FPA. The syntax is: 
. 

dbxenv fpabase [a[0-7]|off] 





If FPA disassembly is disabled (if fpaasm is off) its value is ignored. Other- 
wise, its value is interpreted as follows: 

value in aO . . a7: 

Long move instructions that use the designated address register in 
base+short displacement mode address the FPA, and are disassembled 
using FPA assembler mnemonics. Note that this is independent of the 
actual run-time value of the register. 

value = off: 

All based-mode FPA instructions are disassembled and single-stepped 
as move instructions. 

The default value of fpabase is of f , which designates no FPA base regis- 
ter. 

1. The FPA registers $ f paO . . $f pa31 are recognized and can be used in 
arithmetic expressions or modified in set commands. This extension only 
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B.2. Example of FPA 
Disassembly 



applies to machines with FPA’s. Note that if an FPA register is used in an 
expression or assignment, its type is assumed to be double precision. 

2. FPA registers can be displayed in single precision using the /f display for- 
mat. Double-precision values are displayed using the /F display format. 



Consider the following FORTRAN program: 



/ 


\ 


program example 




print *, f (1 . 0/1.0) 




end 




function f(x,y) 




f = atan (x/y) 




return 




end 




v 


j 



Assume that this program has been compiled with the -g option into the file 
example. On a machine with an FPA, we could disassemble the function f as 
shown below. Note that the FORTRAN intrinsic ATAN is directly supported by 
the FPA instruction set and compiler. 







% dbx a 


. out 


(dbx) stop in f 


(1) stop in f 


(dbx) run 


Running: a. out 


stopped in f at line 5 in file "example . f" 


5 


f = atan (x/y) 


(dbx) &$pc/8i 


f +0x12 


movl a66(0xc),a0 


f +0x16 


fpmoves aO0,fpaO 


f+Oxlc 


movl a60(Ox8),aO 


f +0x20 


fprdivs aO0,fpaO 


f +0x2 6 


fpmoves fpaO, a60 (-Oxc) 


f +0x2e 


fpmoves a60 (-Oxc) , fpal 


f +0x36 


fpatans fpal, fpal 


f +0x40 


fpmoves fpal,a6@ (-0x8) 


v 


j 



FPA disassembly can be disabled by setting the debugger variable f paasm to 
off. This causes dbx to disassemble FPA instructions as long moves to 
addresses on the FPA page: 
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r 

(dbx) dbxenv fpaasm off 
(dbx) &f+0x!2/10i 






f +0x12 


movl 


a60 (Oxc) , aO 




f+0xl6 


movl 


aO0, OxeOOOOcOO : 1 




f +0xlc 


movl 


a60 (0x8) , aO 




f +0x20 


movl 


aO0, 0xe0000600 : 1 




f+0x26 


movl 


Oxe 000 Oe 00:1, a60(-Oxc) 




f+0x2e 


movl 


a60 (-Oxc) , Oxe 0000c 08:1 




f +0x36 


movl 


#0x41, 0xe0000818:l 




f +0x40 


movl 


0xe0000e08 :l,a60(-Ox8) 


J 



When tracing a more complex program, one may occasionally want to step into a 
routine that has been compiled with optimization on. In such routines, the com- 
piled code often addresses the FPA page by using base+short offset addressing. 
Such code can be difficult to recognize unless it is known ahead of time that a 
particular address register is being used to address the FPA. This situation can be 
identified by the presence of an instruction that loads the address of the FPA page 
(OxeOOOOOOO) into an address register before doing any floating-point arithmetic. 



For example, here is a disassembly of the beginning of an optimized FORTRAN 
routine compiled with the -O and -f fpa options: 



r 




(dbx) &ddot_/7i 
ddot : link 


a6, #-0x2a0 


ddot +0x4 : moveml 


#<d2, d3, d4, d5, d6, d7 , a2, a3, a4, a5>, sp0 


ddot +0x8 : lea 


eOOOOOOO : 1, a2 


ddot +0xe : movl 


a20 ( 0xe20 ) ,a60 (-0x278) 


ddot_+0xl4: movl 


a20 ( 0xe24) , a60 (-0x274) 


ddot_+0xla : movl 


a20 (0xe28) ,a60 (-0x270) 


ddot +0x20: movl 


a20 (0xe2c) ,a60 (-0x26c) 


S, 


J 



dbx does not know which register (if any) is being used to address the FPA in a 
given sequence of machine code. However, you may set the dbxenv variable 
f pabase to designate an MC68020 address register as the FPA base register. 
In this example, we note that the compiler has loaded the address of the FPA 
page into register a2, and we then designate a 2 as the FPA base register to 
obtain the following: 



(dbx) dbxenv fpabase a2 
(dbx) &ddot_/7i 

ddot_: link a6,#-0x2a0 

ddot_+0x4 : moveml #<d2 , d3 , d4 , d5, d6, d7 , a2 , a3, a4, a5>, sp@ 

ddot_+0x8 : lea eOOOOOOO : 1, a2 

ddot_+0xe: fpmoved@2 fpa4, a60 (-0x278) 

ddot_+0xla: fpmoved02 fpa5, a60 (-0x270) 

ddot_+0x26: fpmoved02 204ce:l,fpa5 

ddot_+0x36: fpmoved02 204ce:l,fpa4 
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B.3. Examples of FPA FPA data registers can be displayed using a syntax similar to that used for the 

Register Use MC6888 1 coprocessor registers. Note that unlike the MC6888 1 registers, FPA 

registers may contain either single-precision (32-bit) or double-precision (64-bit) 
values; MC6888 1 registers always contain extended-precision (96-bit) values. 



For example, if f paO contains the single-precision value 2.718282, we may 
display it as follows: 



r 






\ 


(dbx) 


&$fpaO/f 






fpa3 


0x402df 855 


+2 . 718282e+00 




< 






J 



Note that the value is displayed in hexadecimal as well as in floating-point nota- 
tion. 



A double-precision value may be displayed using the /F format. For example, if 
fpaO contains the double-precision value 2.718281828, we may display it as 
follows: 



(dbx) 


&$fpaO/F 






fpaO 
< 


0x4005bf Oa 0x8b04919b 


+2.71828182800000e+00 





Note that it is important to use the correct display format; attempting to display a 
double-precision value in single precision will usually produce meaningless 
results. 

FPA registers can also be used in set commands and in arithmetic expressions. 
Since dbx cannot tell whether the value in an FPA register is single- or double- 
precision, dbx provides two sets of names to refer to FPA registers. The names 
$fpaO . . $fpa31 always cause the contents of the register to be interpreted as 
a double-precision value; the names $fpaOs . . $fpa31s cause interpretation 
as a single-precision value. Thus, the commands 



r 







(dbx) 


set $fpaOs = 1.0 




(dbx) 


set $fpa0 =1.0 




V 




j 



store different values in fpaO. 
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FPA Assembler Syntax 



C.l. Instruction Syntax 



This appendix describes the Floating-Point Accelerator (FPA) support extensions 
to as included in Sun software release 3. 1. 

The extensions to as are described in general, with discussions of two-, three-, 
and four-operand instruction examples. Some instructions covered separately 
don’t follow the formats described at the beginning of the appendix. The appen- 
dix includes restrictions and potential errors, followed by a summary of sup- 
ported floating-point instructions. 



The general format for floating-point instructions is 




where 



f p indicates an FPA instruction. 
op is the opcode name. 

t is the operand type, either single (s) or double (d). 

The @A part of the instruction is optional. When present, A specifies the address 
register which contains the base address for the FPA and can be in the range 0..7. 
If this form is used, a previous instruction must load the FPA address 
(OxeOOOOOOO) into the specified address register. 

If @A is not present, then absolute long addressing is used to refer to the FPA. 
This form is more efficient for short routines. 

Depending on the instruction, there may be from zero to four operands specified. 
The operands can be any of the following forms : 

□ Any MC68020 effective address, with the exception that absolute short 
addresses are not allowed for double-precision values. 



□ If either of the data register or the address register is used to hold a double- 

precision value, then the value will be in a register pair and both registers, 
separated by a colon, must be specified in the instruction. For example: 
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C.2. Register Syntax 



C.3. Operand Types 



C.4. Two-Operand 
Instructions 



The only exception to this rule is the fplt od instruction (convert integer 
to double-precision value). 

□ In some instructions (command register type) it is possible to specify that the 
register is in constant RAM. The syntax used for this case is %n, where n 
is a register number in the range 0 to 5 1 1. 

The 32 floating-point data registers are designated f paO , f pal , . . . , 
fpa31. The supported control registers are: 



Hardware 


Software 


MODE3_0 


fpamode 


WSTATUS 


fpastatus 



as supports three floating-point operand types: 

□ s for single-precision operands, 

o d for double-precision operands. 

□ 1 for 32-bit integer operands, used for integer to floating-point conversions. 



Opcodes such as add, subtract, multiply, divide, negate, absolute value, square 
root, conversion from integer to floating-point, conversion from single to double 
(and vice versa) are all represented as: 



r 


*\ 


fp opt X r fpa n 




v 


.j 



where t= s or d, and X is any valid MC68020 effective address for an operand 
or is an FPA data register. 

If X is an FPA register which is in the constant RAM, then it can be in the range 
0 to 51 1. If it is not in constant RAM, then it is one of the 32 FPA data registers. 
When X is an FPA register, then f pare is one of the 32 floating-point data regis- 
ters. If X is an effective address, then f pan is one of the FPA registers in the 
range 0 to 15. The following are examples of such instructions: 



/ 




^ 




Instruction 


Computes 


fpnegs 


<effective addr ess>, fpal 




fpsqrd 


<effective address>, fpa2 




fpsubs 


fpal, fpa2 


fpa2 <— fpa2 - fpal 


fprsubs 


fpal, fpa2 


f^>a2 <- fyal - fjpa2 


fpdivs 


dO, fpa2 


fpa2 <— fpa2 / dO 


fprdivs 


dO, fpa2 


f)>a2 <— dO / fpa2 


V 







In the above examples f pr subs and fprdivs are the reverse subtract and 
reverse divide operators, respectively. 



A 
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The opcodes for sine, cosine, atan, e~x, e~x -1, ln(x), 

In ( 1+x ) and s inco s ( x ) are all supported as command register type instruc- 
tions: 



r 


\ 


fpopt fpax, fpart 






j 



where t = s or d. 

fpaac is either a floating-point register or a register in the constant RAM (which 
is specified as % number). For the sine os instruction, the destination operand 
is actually a register pair: 




C.5. Three-Operand 
Instructions 



where f pac is the cosine’s destination and f pas is the sine’s destination. 

The opcodes +,-,*,/ are supported in extended and command register forms as 

/ \ 

tpop3t X , fpa m, fpa n 

s > 



where t = s or d and X is an <effective address>for an extended instruction or a 
floating-point register for a command register type of instruction. 

In the command register form , X and f para can indicate a register number in the 
constant RAM. That is, they can either be in the range 0 to 5 1 1 or in the range 0 
to 31. In the extended instruction form, fpa m and fpa/a must be in the range 
0 to 15. In the above format the position of X and f pam can be exchanged for 
the commutative operators add and multiply (the result of the operation remains 
the same). 



For example. 




can be represented by either of the following forms: 





fpadd3s 


<effective address>, fpal. 


fpa2 


\ 




fpadd3s 


fpal/ <effective address>. 


fpa2 


J 



The same rule applies to subtract and divide operations. However, they are not 
commutative, so different answers result from each order. For example, 




must be coded as: 



- 








fpsub3s <effective address>, fpal, fpa2 




v 




J 
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whereas 




C.6. Four-Operand 
Instructions 



In the extended and command register formats there are pivot instructions of the 
form: 

\ 

fp opt X, fpax, fpay, fpan 

J 



where fpan is the destination floating-point data register and t = s or d, and X 
is an effective address or a floating-point register. 

In the extended form, the positions of X and f pay can be exchanged for both 
single- and double-precision types of instructions. In single-precision extended 
form, it is possible for two of the four operands to be effective addresses. This is 
in general either the first and third or the second and third operands. 

In the command register form, f pax and f pay can be replaced by %x and %y 
indicating register numbers x and y in the constant RAM. 

For four-operand instructions, f pax, f pay and fpan can each be in the range 
0 to 15, when X is an effective address. If X is an FPA register, then X and 
fpan must be in the range 0 to 3 1 and f pax and f pay can either be in the 
range 0 to 5 1 1 (designating a location in constant RAM) or else in the range 0 to 
31. 

These pivot instructions are rather complicated and will be dealt with com- 
pletely. The following shows the forms of each operation, the assembly code 
equivalent to each form, a generalization of the assembly instruction and the 
sequence of operations equivalent to the pivot instruction. 
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— 


Instruction 


N 

Meaning 


fpma { s , d} 


<effective address>, reg2, reg3, regl 


regl <- reg3 + (reg2 * operand) 


fpma{s, d} 


reg2, reg3 f <effective address>, regl 


regl <— operand + (reg3 * reg2) 


fpma { s , d ) 


reg4, reg2, reg3, regl 


regl <— reg3 + (reg2 * reg4) 


fpmas 


<eal>, reg2, <ea2>, regl 


regl <- operand2 + (reg2 * operandl) 

J 




The fpma instruction, where m 
be generalized as 


stands for multiply, and a stands for add, can 







fpmaf X, fpax, fpay, fpan 

J 



where t is s or d, and X is an effective address>or one of the floating-point 
data registers. In the extended type of instruction, the positions of X and f pay 
can be exchanged. Also, for single precision either the first and third operands or 
the second and third operands can be effective addresses. Note that, for example, 




is equivalent to the following sequence of instructions 



r 






fpmul3s 


dO, fpal r temp 




fpadd3s 


temp f fpa2, temp 




fpmoves 


temp f fpa3 




c 







where temp is a temporary register. 



r 


Instruction 


Meaning 


\ 


fpms { S/ d} 


<effective address> f reg2, reg3 f regl 


regl <— reg3 - (reg2 * operand) 




fpms {s, d] 


reg2, reg3, <effective addres s>, regl 


regl <- operand - (reg3 * reg2) 




fpms { s, d) 


reg4 f reg2, reg3 f regl 


regl <— reg3 - (reg2 * reg4) 




fpmss 


<eal>, reg2 r <ea2> f regl 


regl 4- operand2 - (reg2 * operandl) 


.J 



The f pms instruction, where m stands for multiply, and s stands for subtract, 
can be generalized as 




where t is s or d, and X is an <effective address> or one of the floating-point 
data registers. In the extended type of instruction, the positions of X and fpay 
can be exchanged. Also, in single-precision two-memory instructions, either the 
first and third operands or the second and third operands can be effective 
addresses. Note that, for example. 



r 


A 


fpmss fpal, fpa2, dO , fpa3 
v 


j 
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is equivalent to the following sequence of instructions 



/ 

fpmul3s 


fpal. 


fpa2, temp 


\ 


fpsub3 s 


temp. 


dO, temp 




fpmoves 


temp. 


fpa3 










j 



The fpmr instruction, where m stands for multiply, and r stands for reverse 
subtract, can be generalized as 

. 

fpmrf X r fpax, fpay, fpan 

4 



where t is s or d, and X is an <effective address> or one of the floating-point 
data registers. In the extended type of instruction, the positions of X and fpay 
can be exchanged. 



r 


Instruction 


Meaning 


\ 


fpmr {s,d} 


<effective address>, reg2, reg3, regl 


regl <— (-reg3) + (reg2 * operand) 




fpmr { s , d } 


reg2, reg3, <effective address>, regl 


regl <r- (-operand) + (reg3 * reg2) 




fpmr {s,d} 


reg4, reg2, reg3, regl 


regl <- (-reg3) + (reg2 * reg4) 




fpmrs 

< 


<eal>, reg2, <ea2>, regl 


regl <— (-operand2) + (reg2 * operandl) 


) 



In single-precision extended form either the first and third operands or the second 
and third operands can be effective addresses. Note that, for example, 




is equivalent to the following sequence of instructions: 









fpmul3s 


dO, fpal, temp 




fpsub3s 


fpa2, temp, temp 




fpmoves 


temp, fpa3 




V 







The f pam instruction, where a stands for add, and m stands for multiply, can 
be generalized as 

f v 

fpan \t X r fpax, fpay, fpan 

s / 



where t is s or d, and X is an <effective address> or one of the floating-point 
data registers. In the extended type of instruction, the positions of X and fpay 
can be exchanged. 



#sun 

V microsystems 



Revision A of 12 September 1986 










Appendix C — FPA Assembler Syntax 81 



r 


Instruction 


Meaning 


A 


fpam{s, d} 


<effective address>, reg2, reg3, regl 


regl <— reg3 * (reg2 + operand) 




fpam{ s, d} 


reg2, reg3, <effective address>, regl 


regl <r- operand * (reg3 + reg2) 




fpam{s, d} 


reg4, reg2, reg3, regl 


regl <- reg3 * (reg2 + reg4) 




fpams 

v 


<eal>, reg2, <ea2>, regl 


regl <- operand2 * (reg2 + operand 1) 


V 



In single-precision two-memory instructions, either the first and third operands or 
the second and third operands can be effective addresses. Note that, for example. 



/ 

fpams 

St 


fpal, fpa2, fpa3, fpa4 


_ 

j 


is equivalent to the following sequence of instructions: 


r 

fpadd3s 


fpal, fpa2, temp 




fpmul3s 


temp, fpa3, temp 




fpmoves 

St 


temp, fpa4 


. . j 


The f psm instruction, where s stands for subtract, and 


m stands for multiply, 


can be generalized as 




r 






fpsmf 


X, f par, fpa y, f pa n 








-J 



where t is s or d, and X is an effective address or one of the floating-point data 
registers. In the extended type of instruction, the positions of X and f pay can 
be exchanged. The special cases for single-precision instructions are that either 
the first and third operands or the second and third operands can be effective 
addresses. 



r 




Instruction 


Meaning 





fpsm{ s, d} 


<effective address>, reg2, reg3, regl 


regl <— reg3 * (reg2 - operand) 




fpsm{s,d} 


reg2. 


reg3, <effective address>, regl 


regl <- operand * (reg3 - reg2) 




fpsm{s,d) 


reg4. 


reg2, reg3, regl 


regl <- reg3 * (reg2 - reg4) 




fpsm{s,d} 


reg2. 


<effective addr ess>, reg3, regl 


regl «- reg3 * (-reg2 + operand) 




fpsm{s,d} 


reg2. 


reg4, reg3, regl 


regl <r- reg3 * (-reg2 + reg4) 




fpsms 


<eal>, 


reg2, <ea2>, regl 


regl <— operand2 * (reg2 - operandl) 




fpsms 


reg2, 


<eal>, <ea2>, regl 


regl <— operand2 * (-reg2 + operandl) 


J 



Note that, for example, 



fpsms dO, fpal, fpa2, fpa3 

J 



is equivalent to the following sequence of instructions: 
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— 




N 


fpsub3s 


dO, fpal, temp 




fpmul3s 


temp, fpa2, temp 




fpmove s 


temp, fpa3 











C.7. Other Instructions 



Other special instructions are listed below. In each of them the last operand is 
also the destination, except for tst, cmp and mcmp where fpastatus is 
the implied destination. X is either an effective address or an FPA data register 
and t is either s or d for all instructions except fpmovef, where t can be s, 
d, or 1. 



Table C- 1 Other Instructions 



Mnemonic 


Syntax 


Operation Name 






nop 


fptst / 


X 


operand compare with zero 


fpcmpf 


X, fpam 


register m compare with operand 


fpmcmpf 


X, fpam 


register m compare magnitude with operand 


fpmovef 


fpam, fpan 


move floating-point registers 


fpmove2f 


fpam, fpan 


2x2 matrix move 


fpmove3f 


fpam, fpan 


3x3 matrix move 


fpmove 4 f 


fpam, fpan 


4x4 matrix move 


fpdot2f 


fpax, f pa y, fpan 


fpan <— f par*fpay + 
(fpar+7) * (fpa^+7) 


fpdot3f 


fpax, f pay, fpan 


fpan <— f pajc*fpay + 
(fpa x+1) * (fpay+7) + 
(fpax+2) * (fpa y+2) 


fpdot4f 


fpax, fpa y, fpan 


fpan <— f pajc*fpa> + 

( £paLX+l)*(fp 2 Ly+l ) + (fpax-f 2)*(fpa;y*2) + 
(fpax*f3)*(fpa^^J) 


fptran2/ 


fpam, fpan 


transpose 2x2 matrix 


fptran3f 


fpam, fpan 


transpose 3x3 matrix 


fptran4f 


fpam, fpan 


transpose 4x4 matrix 


fpmove 


fpamode, <ea> 


read mode register 


fpmove 


<ea>, fpamode 


write to mode register 


fpmove 


fpastatus, <ea> 


read status register 


fpmove 


<ea>, fpastatus 


write to status register 


fpmovef 


fpam, <ea> 


read a floating-point data register 


fpmovef 


<ea>, fpan 


write to a floating-point data register 
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C.8. Restrictions and 
Errors 



C.9. Instruction Set 
Summary 



as reports an invalid operand error in double-precision instructions when abso- 
lute short addressing or a single data or address register is used. 

as reports a register out of range error for the dot product and matrix move and 
transpose instructions when the register specified does not fall within the 
specified range. 

For most instructions where one operand is an effective address, the register 
range is 0 to 15. If all operands are FPA registers, then the register range is 0 to 
3 1 . For constant RAM registers, the range is 0 to 5 1 1 . as reports an invalid 
operand error when any of these registers are not within the permitted range. 

In the following table, X is any valid FPA register or MC68020 effective address 
(the form ( xxx ) : w is not allowed for double). In some three- or four-address 
instructions the position of the X and one of the FPA registers can be exchanged 
as shown in the fourth column of the following table. 



Table C-2 Floating-Point Instructions 



Instruction 


Syntax 


Operation 


Alternative 


fpnegs 


X, fpan 


negate single 




fpnegd 


X, fpan 


negate double 




fpabss 


X, fpan 


absolute value single 




fpabsd 


X, fpan 


absolute value double 




fpltos 


X, fpan 


convert integer to single 




fpltod 


X, fpan 


convert integer to double 




fpstol 


X, fpan 


convert single to integer 




fpdtol 


X, fpan 


convert double to integer 




fpstod 


X, fpan 


convert single to double 




fpdtos 


X, fpan 


convert double to single 




fpsqrs 


X, fpan 


square single 




fpsqrd 


X, fpan 


square double 




fpadds 


X, fpan 


add single 




fpadd3s 


X, fpam, fpan 


add single 


fpam, X, fpan 


fpaddd 


X, fpan 


add double 




fpadd3d 


X, fpam, fpan 


add double 


fpam, X, fpan 


fpsubs 


X, fpan 


subtract single 




fpsub3s 


X, fpam, fpan 


subtract single 


fpam, X, fpan 


fprsubs 


<ea>, fpan 


reverse subtract single 




fpsubd 


X, fpan 


subtract double 




fpsub3d 


X, fpam, fpan 


subtract double 


fpam, X, fpan 
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T able C-2 Floating-Point Instructions — Continued 



Instruction 


Syntax 


Operation 


Alternative 


fprsubd 


<ea> , f paw 


reverse subtract double 




fpmuls 


X, fpan 


multiply single 




fpmul3s 


X, fpam, fpan 


multiply single 


fpam, X, fpan 


fpmuld 


X, fpan 


multiply double 




fpmul3d 


X, fpam, fpan 


multiply double 


fpam, X, fpan 


fpdivs 


X, fpan 


divide single 




f pdi v3 s 


X, fpam, fpan 


divide single 


fpam, X, fpan 


fprdivs 


<ea> , fpan 


reverse divide single 




fpdivd 


X, fpan 


divide double 




fpdiv3d 


X, fpam, fpan 


divide double 


fpam, X, fpan 


fprdivd 


<ea> , fpan 


reverse divide double 




fpnop 




nop 




fptsts 


X 


single compare with 0 




fptstd 


X 


double compare with 0 




fpcmps 


X, fpam 


single compare 




fpcmpd 


X, fpam 


double compare 




fpmcmps 


X, fpam 


single magnitude compare 




fpmcmpd 


X, fpam 


double magnitude compare 




fpsins 


fpax, fpan 


sine single 




fpsind 


fpax, fpan 


sine double 




fpcoss 


fpax, fpan 


cosine single 




fpcosd 


fpax, fpan 


cosine double 




fpatans 


fpax, fpan 


atan single 




fpatand 


fpax, fpan 


atan double 




fpetoxs 


fpax, fpan 


e~x single 




fpetoxd 


fpax, fpan 


e~x double 




fpetoxmls 


fpax, fpan 


e~x-l single 




fpetoxmld 


fpax, fpan 


e~x-l double 




fplogns 


fpax, fpan 


In (x) single 




fplognd 


fpax, fpan 


In (x) double 




fplognpls 


fpax, fpan 


ln(l+x) single 




fplognpld 


fpax, fpan 


ln(l+x) double 




fpsincoss 


fpax, f pac : fpa s 


fpac <— cosine (x), fpas <— sine(x) 




fpsincosd 


fpax, f pac : f paj 


fpac <— cosine (x) , fpas f— sine (x) 




fpmas 


X, fpax, fpa}, fpan 


fpan (fpax * X) + fpa} 










fpax, X, fpa}, fpan 








fpa}, fpax, X, fpan 








X, fpax, X, fpan 








fpax, X, X, fpan 


fpmad 


X, fpax, fpa}, fpan 


fpan (fpax * X) + fpa} 
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Table C-2 Floating-Point Instructions — Continued 



Instruction 


Syntax 




Operation 


Alternative 










fpax, 


X, fpay, fpan 










fpay. 


fpax, X, fpan 


fpmss 


X, fpax, fpay. 


fpan 


fpan 4 — fpay - (fpax: * X) 


fpax. 


X, fpay, fpan 










fpay, 


fpax, X, fpan 










X, fpax, X, fpan 










fpax, 


X, X, fpan 


fpmsd 


X, fpax, fpay, 


fpan 


fpan 4 - fpay - (fpax: * X) 


fpax, 


X, fpay, fpan 










fpay, 


fpax, X, fpan 


fpmrs 


X, fpax, fpay. 


fpan 


fpan 4 - (fpax: * X) - fpay 


fpax. 


X, fpay, fpan 










fpay, 


fpax, X, fpan 










X, fpax, X, fpan 










fpax. 


X, X, fpan 


fpmrd 


X, fpax, fpay. 


fpan 


fpan 4 — (fpax * X) - fpay 


fpax, 


X, fpay, fpan 










fpay, 


fpax, X, fpan 


fpams 


X, fpax, fpay, 


fpan 


fpan <— (fpax + X) * fpay 


fpax. 


X, fpay, fpan 










fpay, 


fpax, X, fpan 










X, fpax, X, fpan 










fpax, 


X , X , fpan 


fpamd 


X, fpax, fpay, 


fpan 


fpan 4 - (fpax + X) * fpay 


fpax. 


X, fpay, fpan 










fpay. 


fpax, X, fpan 


fpsms 


X, fpax, fpay. 


fpan 


fpan 4 — (fpax - X) * fpay 


fpax. 


X, fpay, fpan 










fpay. 


fpax, X, fpan 










X, fpax, X, fpan 










fpax. 


X, X, fpan 


fpsmd 


X, fpax, fpay. 


fpan 


fpan 4 - (fpax - X) * fpay 


fpax. 


X, fpay, fpan 










fpay. 


fpax, X, fpan 


fpmoves 


<ea> , f pan 




write to a register, single 






fpmoved 


<ea>, fpan 




write to a register, double 






fpmovel 


<ea> , fpan 




write to register, integer 






f pmove s 


fpa/n, <ea> 




read a register, single 






fpmoved 


fpam, <ea> 




read a register, double 






f pmove 2s 


fpa/n, fpan 




2x2 matrix move, single 






f pmove 2d 


fpa/n, fpan 




2x2 matrix move, double 






f pmove 3s 


fpam, fpan 




3x3 matrix move, single 






f pmove 3d 


fpam, fpan 




3x3 matrix move, double 






f pmove 4s 


fpam, fpan 




4x4 matrix move, single 






fpmove4d 


fpam, fpan 




4x4 matrix move, double 
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Table C-2 Floating-Point Instructions — Continued 



Instruction 


Syntax 


Operation 


Alternative 


fpdot2s 


fpax, fpa}?, fpan 


fpan <- fpax* fpa}? + (fpax-fi) * (fpa y+1) 




fpdot2d 


fpax, fpa}?, fpan 


fpan <— fpax* fpa}? + (fpa x+1) * (fpa y+1) 




fpdot3s 


fpax, fpa}?, fpan 


fpan <— fpax* fpa}? + (fpax+7) * (fpa y+1) + 








(fpax+2) * (fpa y+2) 




fpdot3d 


fpax, fpa y, fpa n 


fpan <— fpax* fpa y + (fpax+7) * (fpa y+1) + 








(fpax+2) * (fpa y+2) 




f pdot 4 s 


fpax, fpa y, fpa n 


fpan fpax*fpa}? + (fpax+1) * (fpay-hl) + 








(fpax-f2) * (£pay+2) + (fpax-f3) * (fpa}?-f3) 




fpdot4d 


fpax, fpa}?, fpan 


fpan <— fpax*fpa}? +■ (fpax-f 7) * (fpa}? -f 7 ) + 








(fpaz+2) * (fpa>+2) + (fpax+3) * (fpay+J) 




fptran2s 


fpam, fpan 


transpose 2x2 matrix, single 




fptran2d 


fpa m, f paw 


transpose 2x2 matrix, double 




fptran3s 


fpam, fpan 


transpose 3x3 matrix, single 




fptran3d 


fpam, f pan 


transpose 3x3 matrix, double 




fptran4s 


fpam, fpan 


transpose 4x4 matrix, single 




fptran4d 


fpam, fpan 


transpose 4x4 matrix, double 




fpmove 


fpamode, <ea> 


read the mode register 




fpmove 


<ea>, fpamode 


write on mode register 




fpmove 


fpastatus, <ea> 


read the status register 




fpmove 


<ea>, fpastatus 


write to status register 
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IEEE Appendix Functions 



The IEEE Standard for Binary Floating-Point Arithmetic includes an appendix of 
useful functions which are not usually available in higher-level languages. The 
IEEE test vectors measure compliance with the specification of some of those 
functions. The following paragraphs give examples of double-precision FOR- 
TRAN and assembly language routines that implement those functions in a way 
that passes the IEEE test vectors under conditions described under "Conformance 
Benchmarks" in Chapter 4. The assembly language examples are coded to use 
switched floating-point. These functions, especially the assembly language ones 
that call V... and F... entry points in libc, may not work in future software 
releases. Efficiency was not a coding consideration. 



Copysign: 



— 






— 


.globl 


_copyd_ 


| Works with any kind of floating point! 




_copyd_: 








movel 


sp@ (4) , aO 






mo vein 1 


a0@, dO/dl 


I dO := x. 




movel 


sp0 (8) ,a0 


| Address of y. 




tstb 


aO0 






bmis 


if 






bclr 


#31, dO 






bras 


2f 






1: 








bset 


#31, dO 






2: 








rts 








v 






J 



Logb: 



r 






.globl 


_logbd_ 




_logbd_: 






movel 


sp@ (4) ,a0 




movel 


a0@, dO 




jsr 


Fexpod 




cmpw 


#-0x3ff,d0 




beqs 


3f 




cmpw 


#0x400, dO 




beqs 


2f 




^ 




. > 
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r 1 : 


jsr 


Vf ltd 








“ \ 




rts 












2: 


moveml 


aO0,dO/dl 












bclr 


#31, dO 












lea 


dzero, a0 












jsr 


Vaddd 


I Add 


zero to 


catch signalling NaN. 






rts 












3: 


moveml 


aO0,dO/dl 












bclr 


#31, dO 












tstl 


dO 












bnes 


5f 












tstl 


dl 












beqs 


4f 










5: 


movel 


#-0x3f f , dO 












bras 


lb 










4: 


moveml 


dmone , dO /dl 












lea 


dzero, aO 












jsr 


Vdivd 


1 Set 


up -1/0 


to generate divide by zero signal. 






rts 












dzero : 


. double 


OrO 










dmone : 


. double 


Or-1 










S 












/ 



Scalb: 



r 




\ 


real*8 


function scalbd ( x, y) 




real*8 


x, y, scalid 




scalbd 


= scalid ( x, int(y)) 




end 






.globl 


_scalid_ 




_scalid_: 






movel 


sp0 (4) , aO 




moveml 


aO0,dO/dl 


1 x. 


movel 


sp0 (8) , aO 


| Address of y. 


jsr 


Vscaleid 




rts 










J 
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Fraction part: The fraction part or significant! function, not in the Standard, is 
nonetheless tested by an IEEE test vector. 

real*8 function signifd ( x ) 
real*8 x, logbd, scalbd, s 
s = scalbd (x, -logbd (x) ) 

if ((abs(s) .gt. 0 . OdO ) .and. (abs (s) .It. l.OdO)) then 
1 s = s + s 

if (abs (s) .It. l.OdO) goto 1 

endif 

signifd = s 
end 



Nextafter: 



real*8 function nextd ( x, y) 
real* 8 x, y, t, minnorm, maxnorm 
integer i(2) 
equivalence (i,t) 

parameter ( minnorm = 2 . 225073858507201383d-308 ) 
parameter ( maxnorm = 1 . 7 97 6931348623157d+308) 

if ( (y .ne. y) .or. (x .ne. x) ) then 

handle nans; convert signalling nans 
nextd = x + y 
return 
end if 

t - x 

if ( x .It. y) then 

if (x .It. 0) then 

i (2) = i (2) - 1 

if (i (2) .eq. -1) i(l) = i (1) - 1 
else if (x .gt. 0) then 
i (2) = i (2 ) + 1 

if (i (2) .eq. 0) i (1) = i(l) + 1 
else if (x .eq. 0) then 
i (2) = 1 
i(l) = 0 



endif 

else if ( x 
if 
i (2 ) 



.gt. y) then 
(x .It. 0) then 

- i (2 ) + 1 
if (i (2) .eq. 0) i(l) = i(l) + 1 
else if (x .gt. 0) then 
i (2 ) - i (2) - 1 

if (i (2) .eq. -1) i (1) = i(l) - 3 
else if (x .eq. 0) then 

i (2) = 1 
i(l) « 0 

t = -t 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



c 



end if 
end if 

It would have been preferable to use the following statements to 
generate underflow and overflow signals, but optimization with -0 
eliminates the statements as dead code, 
if (x .ne. t) then 

if (abs(t) .It. minnorm) then 
signal underflow without changing t 
z = minnorm * 0 . ldO 
else if (abs(t) .gt. maxnorm) then 
signal overflow without changing t 

z = maxnorm * 2 . OdO 

end if 

The following code accomplishes the same thing with more trouble, 
but fails with -f68881 -O because of extended register allocation. 
if (x .ne. t) then 

if (abs(t) .It. minnorm) then 
signal underflow without changing t 

t = copyd(min (abs (t) , 4. OdO * (minnorm/3 . OdO) ) , t) 
else if (abs(t) .gt. maxnorm) then 
signal overflow without changing t 

t = copyd (max (abs (t) , maxnorm + maxnorm), t) 

end if 

if (x .ne. t) then 

if (abs(t) .It. minnorm) then 
signal underflow without changing t 

call dummyd (minnorm * minnorm) 
else if (abs (t) .gt. maxnorm) then 
signal overflow without changing t 

call dummyd (maxnorm * maxnorm) 

end if 
end if 
nextd ~ t 
end 

subroutine dummyd ( x ) 

real*8 x 

end 
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SPICE Input Files 95 
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SPICE Input Files 



E.l. tdo 



The following input file models a tunnel diode oscillator; the transient analysis 
output calculation is extremely ill-conditioned, as discussed under "Different 
Numerical Results" in Chapter 2: 



TDO - TUNNEL DIODE OSCILLATOR 
.WIDTH IN=72 
VBIAS 0 2 -120MV 

LS 21 2.5UH 

CS 10 100PF 

GTD 1 0 POLY ( 1 ) 1 0 

+ -3.95510115972848E-17 +1 . 80727308405845E-01 -2 . 93646217292003E+00 

+ +4.12669748472374E+01 -6 . 0 964951 6869413E+02 +6 . 08207899870511E+03 

+ -3. 7345933647 87 68E+04 +1 . 44146702315112E+05 -3 -53021176453665E+05 

+ +5. 340934360847 62E+05 -4 . 56234076434067E+05 +1 . 68527934888894E+05 

-DC VBIAS 0 -600MV -5MV 
-PLOT DC I (VBIAS) (0,5MA) 

-TRAN 5NS 500NS 0 5NS 
-PLOT TRAN V(l) 

-OPT ACCT LIST NODE LVLCOD=2 
-END 
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E.2. mo samp 2 



The following input file models an MOS amplifier, and is used to obtain the 
benchmark results discussed under "Performance Benchmarks" in Chapter 4: 



mosamp2 - mos amplifier - transient 
.WIDTH OUT=80 



.OPTIONS 


ACCT ABSTOL=10N 


VNTOL=10N 


.TRAN 


o.: 


1US 


10US 






Ml 


15 


15 


1 


32 


M 


W=88.9U 


L=25.4U 


M2 


1 


1 


2 


32 


M 


W=12 . 7U 


L=266.7U 


M3 


2 


2 


30 


32 


M 


W=88.9U 


L=25 . 4U 


M4 


15 


5 


4 


32 


M 


W=12.7U 


L=106.7U 


M5 


4 


4 


30 


32 


M 


W=88.9U 


L=12.7U 


M6 


15 


15 


5 


32 


M 


W=44.5U 


L=25.4U 


M7 


5 


20 


8 


32 


M 


W=482.6U 


L=12 . 7U 


M8 


8 


2 


30 


32 


M 


W=88 . 9U 


L=25.4U 


M9 


15 


15 


6 


32 


M 


W=44 . 5U 


L=25.4U 


M10 


6 


21 


8 


32 


M 


W=482. 6U 


L=12 . 7U 


Mil 


15 


6 


7 


32 


M 


W=12 . 7U 


L=106.7U 


M12 


7 


4 


30 


32 


M 


W=88.9U 


L=12 . 7U 


Ml 3 


15 


10 


9 


32 


M 


W=139.7U 


L=12 . 7U 


Ml 4 


9 


11 


30 


32 


M 


W=139.7U 


L=12.7U 


Ml 5 


15 


15 


12 


32 


M 


W=12.7U 


L=207.8U 


Ml 6 


12 


12 


11 


32 


M 


W=54.1U 


L=12 . 7U 


Ml 7 


11 


11 


30 


32 


M 


W=54 . 1U 


L=12 . 7U 


Ml 8 


15 


15 


10 


32 


M 


W=12 . 7U 


L=45 . 2U 


Ml 9 


10 


12 


13 


32 


M 


W=270.5U 


L=12 . 7U 


M20 


13 


7 


30 


32 


M 


W=270.5U 


L=12 . 7U 


M21 


15 


10 


14 


32 


M 


W=254U 


L=12 . 7U 


M2 2 


14 


11 


30 


32 


M 


W=241.3U 


L=12 . 7U 


M2 3 


15 


20 


16 


32 


M 


W=19U 


L=38 . 1U 


M2 4 


16 


14 


30 


32 


M 


W=406.4U 


L=12 . 7U 


M2 5 


15 


15 


20 


32 


M 


W=3 8 . 1U 


L=42.7U 


M2 6 


20 


16 


30 


32 


M 


W=381U 


L=2 5 . 4U 


M2 7 


20 


15 


66 


32 


M 


W=22 . 9U 


L=7 . 6U 


CC 


7 9 


40PF 










CL 


66 i 


0 70PF 











VIN 21 0 PULSE (0 5 INS INS INS 5US 10US) 

VCCP 15 0 DC +15 
VDDN 30 0 DC -15 
VB 32 0 DC -20 

.MODEL M NMOS (NSUB=2 . 2E15 UO=575 UCRIT=49K UEXP=0 . 1 TOX=0.11U XJ=2 . 95U 
+ LEVEL=2 CGSO=1.5N CGDO=l . 5N CBD=4.5F CBS=4.5F LD=2.4485U NSS=3.2E10 
+ KP=2E-5 PHI=0 . 6 ) 

.PRINT TRAN V(20) V(66) 

-PLOT TRAN V (20) V(66) 

-END 
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fmove[sd] fpm,<ea> 99 

f mo vex fp/n,<ea> 99 

flognt, flog2t, floglOt 100 

fsqrtx fpm, fp« (mon) 100 

fsqrtp <ea>,fpn 100 

Binary-to-Decimal Conversion 100 

Decimal-to-Binary Conversion 101 




MC 6888 1 Mask Differences 



fmove [ sd] 



fp m, <ea> 



fmovex fpm,<ea> 



This appendix describes differences between the MC6888 l’s A79J and later 
A93N masks. The A93N mask reflects the MC6888 1 described in Motorola’s 
definitive MC68881 Floating-Point Coprocessor User’ s Manual. 

Youcanrun mc68881version ( 8 ) on a Sun-3 to determine the installed 
MC68881’s mask and approximate clock rate. 

Some of the problems listed below are worked around in software in Sun’s 
Release 3.0 and later, sometimes at a cost in accuracy and performance, relative 
to A93N hardware and software without workarounds. Since one of the problems 
described below has no workaround, and future releases of Sun software may 
have other workarounds removed to improve accuracy and performance, it is a 
good idea for you to upgrade to the A93N-version MC6888 1 when it is available. 

fmove [ sd] f pm, <ea> to single- or double-precision destinations produces 
an incorrect result when the result underflows but rounds back up to the smallest 
normalized number. The result returned is zero instead of the smallest normal- 
ized number. 

This problem might occasionally confound programs that depend on monotoni- 
city: there are single- and double-precision x > 0 such that the computed result 
of x / 2 is 0 but the computed result of x / 4 is > 0. There is no software wor- 
karound for this problem, which may, but probably does not, adversely affect any 
program in which underflow occurs. Underflow can be detected by enabling the 
underflow bit in the MC6888 1 Exception Enable Byte. 

fmovex f pm, <ea> to an extended-precision destination produces an 
incorrect result when the source operand is an extended-precision denormalized 
number. The result returned is zero, an incorrect denormalized number, or 
an incorrect, tiny normalized number. 

The problem can be demonstrated by constructing an appropriate sequence of 
operations on double-precision operands but it is unlikely to be encountered 
unintentionally. There is no software workaround for this problem, but it is 
unlikely to affect any Sun-supplied software since the extended data type is not 
direcdy supported by Sun’s compilers. 
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flognt, flog2t, 
floglOt 



fsgrtx fpm, fpn ( mon ) 



fsqrtp <ea>, fpn 



Binary -to-Decimal Conversion 



The flognt, flog2t, and floglOt instruction produces an excessive 
amount of error in the results for operands lie. As e approaches 0, the results 
become totally dominated by error. 



Sun’s software releases 3.0 and later contain workarounds that avoid this prob- 
lem for programs written in higher-level languages, as the assembly language 
entry Mlogd illustrates: 



f 




\ 


fmoved sp0,fpO 


1 sp0 contains 


the argument x. 


fcmps #0r0.5,fp0 


f jule If 


| Branch if x<= 


=0.5 or x is Nan. i 


fsubl #l,fp0 


flognplx fpO,fpO 


IThis is more 


accurate for x>0.5. 


bras 2f 

*1 _ 


JL • 

flognx fp0,fp0 

9 . 


| This is more 


accurate for x<0.5. 


fmoved fpO,sp0 


| sp0 receives 


the computed logn. 


V 




J 



fsqrtx fpm, fpn (mon) can produce an incorrect result; the failure is 
dependent upon the contents of fpn before the fsqrtx executes. There- 
fore, this problem can be avoided in software by preceding the fsqrtx with an 
instruction which loads fpn with an innocuous value such as 2.0: precede the 
fsqrtx with fmoveb #2 , fpn or an equivalent instruction. 

fsqrtp <ea>, fpn may cause the same problem, and may be avoided in the 
same way. Note that this problem can’t occur for either register-to-register 
fsqrtx f pm, f pm or memory-to-register f sqrt [sdxbwl] <ea>,fpn. 

Sun’s software releases 3.0 and later contain workarounds that avoid this prob- 
lem for programs written in higher-level languages. 

a. When die result is an exact power of 10, the packed BCD mantissa is 
$C.OOO....OOO instead of $1.000.. .000 (note non-BCD digit). 

b. The result will have the incorrect exponent sign when the decimal 
exponent is greater than +999. 

c. When the "K" parameter is in the range 0 thru 17, die OPERR (IOP) 
exception is incorrectly signalled. 

d. When the magnitude of the decimal exponent is greater than 999, the 
OVFL exception is incorrectly signalled (OPERR should be signalled). 

e. When the source operand is an extended-precision denormalized 
number, the numeric value of the result is correct but the decimal round- 
ing boundary is incorrect 

f. The inexact exception is set even when conversion is exact 

g. Conversions with a "K" parameter of 0 do not function as specified, 
the result is the same as K=+l. 
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Sun’s releases 3.0 and later use the MC6888 l’s binary -to-ASCII conversion 
instructions only in the adb and dbx debuggers. These problems are worked 
around there. 

Decimal-to-Binary Conversion f opp <ea>, f pm of some exact powers of ten contains an error of 1 bit. 

These instructions are not used in Sun’s releases 3.0 and later. 
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G.l. Introduction 



Language-Specific Constructs 



Access to Special Instructions 



Register Allocation 




Assembly-Level In-line Expansion 



A simple in-line expansion facility lets you integrate assembly routines into 
higher-level routines written in C, Pascal, or FORTRAN. 

The peephole optimizer c2 has been modified to optimize code sequences con- 
taining inline-expanded routines. In many cases, the resulting code is compar- 
able in quality to that compiled for constructs directly supported in the common 
code generator. 

In-line expansion of procedure calls is an important optimization strategy, pri- 
marily because it exposes opportunities for other optimizations. It is also useful 
as a "software Swiss Army knife" in applications that depend on routines written 
in assembly language, either for performance reasons or for access to machine 
instructions not otherwise available in high-level languages. Several such appli- 
cations are described in the following sections. 

The constructs of a given language frequently cannot be supported in a common 
backend with a reasonably small operator set in multilingual environments based 
on a common code generator. For example, Pascal sets and strings, and FOR- 
TRAN complex numbers are not naturally supported in a compiler system based 
on the Portable C Compiler (pcc), which includes none of these types. Attempts 
to support all such language-specific constructs in a single code generator typi- 
cally result in large, unwieldy code generators containing large amounts of code 
intended only for a specific language. In-line expansion of language-specific 
library routines presents an attractive alternative. 

A programmer frequently requires access in system software to instructions not 
normally produced by the code generator. The usual technique is to call a hand- 
coded assembly routine which uses the desired instruction. Unfortunately, the 
procedure call and return overhead may be unacceptable if this occurs in a region 
of high execution frequency. 

As noted earlier, in-line expansion is an important optimization strategy, pri- 
marily becauses it exposes opportunities for improved register allocation and 
other optimizations. This is important in the Sun Workstation environment, for 
several reasons. 

First, the Sun MC68000 family calling sequence requires that all argument 
values be passed on the stack. This incurs memory traffic that becomes unneces- 
sary when the body of a called procedure is integrated into the caller’s text 
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G.2. User Interface 



Implementation 



Second, the calling sequence requires that all function results be returned in the 
main processor registers dO and dl. This is particularly annoying for functions 
returning floating-point values, since the calling routine cannot normally assume 
that the called routine uses a specific floating-point processor. Consequently, the 
calling sequence typically incurs a significant amount of traffic between the main 
processor registers and the floating-point processor registers. For simple opera- 
tions, the cost of moving data in and out of external registers can easily exceed 
the cost of the operation itself. Experience shows that such operations are used 
frequently. 

In the C, FORTRAN, and Pascal compilers, filenames ending with the suffix 
" . il" are assumed to contain inline-expandable assembly routines. In addition, 
the Pascal compiler uses the same mechanism to expand standard procedures and 
functions. 

In-line expansion is divided into two tasks. The expansion itself is a straightfor- 
ward text substitution, with little knowledge of the details of the assembly 
language. Most of the work of reducing calling sequence overhead is done by 
the peephole optimizer c2. This partitioning of tasks is convenient, since most 
of the information required is already present in the program representation built 
by c2 for other optimizations. 

For illustrative purposes, consider the FORTRAN function cexp() or complex 
exponential. Assuming the real functions sin(), cos(), and exp() as primi- 
tives, the cexpQ function may be implemented by the following C routine: 
. 

typedef struct { 

double real,imag; 

} dcomplex; 

void c_exp ( r , z ) 

dcomplex *r, *z; 

{ 

register double expx; 
double exp(), cos(), sin(); 
expx = exp(z->real) ; 
r->real = expx * cos (z->imag) ; 
r->imag = expx * sin (z->imag) ; 

} 

^ > 



Straightforward compilation with MC6888 1 code generation (- f 6 8 8 8 1) 

enabled produces the following translation: 
. 



• text 

I #PROC# 04 


.globl 


_c_exp 


_c_exp : 


link 


a6, #0 


addl 


#-LF12,sp 


moveml 


#LS12, sp0 


fmovem 


#LSS12,a6@ (-LFF12) 
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movl 


a60 (Oxc) , aO 


movl 


aO0 (0x4) , sp0- 


movl 


aO0, sp0- 


jbsr 


_exp 


addqw 


#0x8, sp 


movl 


dl, sp0- 


movl 


d0, sp0- 


fmoved 


sp0 + , fp7 


movl 


a60 (Oxc) , aO 


movl 


aO0 (Oxc) , sp0- 


movl 


aO0 (0x8) , sp0- 


jbsr 


_cos 


addqw 


#0x8, sp 


movl 


dl, sp0- 


movl 


dO, sp@- 


fmoved 


Sp0+, fpO 


fmulx 


fp7,fp0 


movl 


a60 (0x8) , aO 


fmoved 


fpO , aO0 


movl 


a60 ( Oxc) , aO 


movl 


aO0 (Oxc) , sp0- 


movl 


aO0 (0x8) , sp0- 


jbsr 


__sin 


addqw 


#0x8, sp 


movl 


dl, sp0- 


movl 


dO, sp0- 


fmoved 


sp0 + , fpO 


fmulx 


fp7 , fpO 


movl 


a60 (0x8) , aO 


fmoved 


f p0 , aO0 (0x8) 


fmovem 


a60 (-Oxc) ,#0xl 


unlk 

rts 


a6 


LF12 = 


0 


LS12 = 


0x0 


LFF12 = 


* 12 


LSS12 = 


= 0x1 


LP12 = 
.data 


0x10 



Note that most of this code is occupied with passing arguments on the stack, and 
moving function results from dO /dl to f pO. Even the latter involves stack 
traffic, since the MC6888 1 does not support direct moves between floating-point 
registers and register pairs on the main processor. Additional overhead (not 
shown) exists for similar reasons in the called library routines. In-line expansion 
can do much to alleviate such problems, as will be shown in the following sec- 
tions. 
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G.3. In-line Expansion Pass 



Inline, die program that actually does in-line expansion, is little more than a 
glorified sed script It knows nothing about the syntax or semantics of the target 
machine assembly language, other than the form of a call instruction. Inline 

is invoked by the cc, pc, and f 77 compilers as 
- 

/usr/ lib/inline [sourcefile] [-o outputfile ] [-i inlinefile ] 

V ) 



Inline expands call instructions in sourcefile using routine definitions from 
one or more inlinefile^. If no sourcefile is specified, the default source is 
stdin. If no outputfile is specified, the default output is stdout. If no 
inlinefile is specified, inline behaves essentially like /bird cat. Note that 
nested expansions are not supported, but may be implemented using pipes. 

Each inlinefile contains one or more labeled assembly language routines of the 
form: 

. 

.inline name,argsize 

instructions 

• end 

^ 



where the instructions constitute an in-line expansion of the named routine, and 
argsize is the number of bytes of arguments expected. The routine must observe 
the following restrictions: 

[1] Registers a0-al/d0-dl/fp0-fpl/fpa0-fpa3 may be used freely. 

[2] Other registers must be saved on entry and restored on exit 

[3] Results are returned in dO or dO/dl. 

[4] Arguments must be explicitly deleted from the stack. In general, this should 
be done using autoincrement addressing. 

The optimizations performed in c2 assume that in-line routines are coded in a 
way that makes the lifetimes of stack temporaries explicit On the MC68000 
family of processors, this can be done by using autoincrement addressing to pop 
incoming arguments from the stack. For the preceding example, on the 
MC6888 1 you could code double-precision versions of sin(), cos(), and 
exp() as follows: 



r 






N 




.inline 


_cos, 8 






fcosd 


sp0+, fpO 






fmoved 


fpO, sp@- 






movl 


sp0+, dO 






movl 


sp0+, dl 






.end 








. inline 


_sin, 8 






f sind 


sp0+, fpO 






fmoved 


fpO r sp0- 




\ 






J 
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movl 


sp0+,dO > 


movl 


sp0+,dl 


.end 




. inline 


__exp, 8 


fetoxd 


sp0+, fpO 


fmoved 


fpO, sp0- 


movl 


sp0+,dO 


movl 


sp0+,dl 


.end 


/ 



When inline is run on the preceding compiled code using these routines, the 
result is: 



.text 

.globl _c_exp 

c_exp : 

link a6,#0 

addl #-LFl2,sp 

moveml #LSl2,sp0 

fmovem #LSSl2,a60 (-LFF12) 

movl a60(Oxc),aO 

movl a0@ (0x4) , sp0- 

movl aO0,sp0- 

fetoxd sp0+,fpO 

fmoved fpO,sp0- 

movl sp0+,dO 

movl sp0+ f dl 

subql #8,sp 

addqw #0x8 , sp 

movl dl,sp0- 

movl dO,sp0- 

fmoved sp0+, fp7 

movl a60(Oxc) f aO 

movl aO0 (Oxc) , sp0- 

movl aO0 (0x8) , sp0- 

fcosd sp0+,fpO 

fmoved f p0 , sp0- 

movl sp0+,dO 

movl sp0+ f dl 

subql #8,sp 

addqw #0x8, sp 

movl dl,sp0- 

movl dO,sp0- 

fmoved sp0+,fpO 

fmulx fp7,fp0 

movl a60 (0x8) , aO 

fmoved fpO,aO0 

movl a60 (Oxc) , aO 

movl aO0 (Oxc) , sp0- 

movl aO0 (0x8) , sp0- 

fsind sp0+,fpO 
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fmoved 


fpO, sp0- 


movl 


sp0+,dO 


movl 


sp0+, dl 


subql 


#8, sp 


addqw 


#0x8, sp 


movl 


dl, sp0- 


movl 


d0, sp0- 


fmoved 


sp0+, fpO 


fmulx 


fp7 , fpO 


movl 


a60 (0x8) , aO 


fmoved 


fpO, aO0 (0x8) 


fmovem 


a60 (-Oxc) ,#0xl 


unlk 


a6 


rts 





G.4. Peephole 

Optimizations 



This code looks terrible, but it has several desirable attributes. 

1. It will execute correctly in its present form, and may even be slightly faster 
than the original code, although it is considerably larger. 

2. All stack traffic is explicit, rather than being hidden in the semantics of the 
procedure call. 

3. The patterns of stack overhead are regular and may be readily optimized by 
local transformations. In particular, no elaborate analysis is necessary to 
determine lifetimes of stack temporaries. 



The symmetric use of autoincrement and autodecrement addressing is taken by 
c2 to indicate the lifetime of the stack copy of an argument. In particular, if 
there are no other uses of the stack copy, the copy need not be created, and the 
argument may instead be loaded into a register or used directly in the source field 
of an instruction in the expanded routine. This is typical of a class of optimiza- 
tions in c2 which are similar to copy propagation, and generally attempt to 
reduce the height of the runtime stack. These optimizations depend on a program 
representation in which uses, definitions, and lifetimes of registers are explicit. 
For example, consider the following code pattern: 



— 
(1) move x,sp0- 




(2) <op> sp0+, . . . 


> 


If possible, this pattern is rewritten as 


(1) (deleted) 


A 


(2) <op> x , . . . 




J 



This optimization is feasible only if x is not modified between (1) and (2). Other 
requirements include: side effects of (1) (including condition codes) must not be 
used, the value on the top of the stack must not be used, and the stack pointer 
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itself must not be used or set If any of these restrictions are not met, the code is 
left unchanged. 



Another similar pattern is the following: 


(1) move Jt, sp0- 


N 


(2) move sp0+ ,rn 

< 


J 


If possible, this pattern is rewritten as 


( 1 ) move x t rn 




(2) (deleted) 

V 





This optimization requires that the register rn be neither used nor set in the 
interim. Care must also be taken to leave the code unchanged if the condition 
codes are live after either instruction (1) or instruction (2). 

Applying these and other c2 optimizations to the example of the previous sec- 
tion yields the following: 

V 



.text 

.globl 


_c_exp 


link 


a6,#-12 


fmovex 


f P 7, a60 (-12) 


movl 


a6@ (12), aO 


fetoxd 


aOG , fp7 


fcosd 


a0@ (8) ,fpO 


fmulx 


fp7 , fpO 


movl 


a60 (8) , aO 


fmoved 


fpO / aO0 


movl 


a60 (12) ,a0 


f sind 


a0@ (8) , fpO 


fmulx 


fp7 , fpO 


movl 


a6G (8) ,a0 


fmoved 


fpO, a06 (8) 


fmovex 


a60 (-12) , fp7 


unlk 


a6 


rts 





Although this code is not optimal, the improvement over the initial version is 
substantial. In addition, the final code is smaller than that which would be 
obtained by optimizing the original unexpanded code. Thus a gain in speed is 
obtained without incurring an increase in code size. 
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G.5. Using Sun’s 

Predefined . il Files 



Faster Execution 



Smaller Executable Files 



Starting with Release 3.1, Sun provides five inline expansion template files: 
/usr/ lib/ f switch . il 
/ usr/ lib/ f soft . il 
/ usr/lib/ f sky . il 
/usr/lib/f 68881.il 
/usr/ lib/ f fpa . il 

These files may be optionally used to replace calls to procedures in libF7 7 . a 
and libm . a with code that will either execute faster or provide smaller execut- 
able files. 

The procedure is simple; replace a command like 

"'i 

ill -0 -ffpa -c cx.f 





■with 

— 
ill -0 -ffpa -c cx.f /usr/lib/ffpa . il 
^ > 



Then any code generated to call routines redefined in / usr/ lib/f fpa . il 
will be appropriately expanded inline. More than one . il file can be specified; 
in that case, the first definition encountered of an inline-expanded procedure will 
supersede any following definitions of the same procedure. 

None of Sun’s compilers use the f * . il files from /usr/lib except when 
specifically requested, pc always automatically uses another file pc2 . il 
which doesn’t affect floating-point subroutines. 

The /usr/lib/f* . il files’ primary application is to accelerate calculations 
involving the complex and doublecomplex data types in FORTRAN. For 
-ffpa, -f 68881, and -f sky, intensive complex arithmetic may be twice as 
fast with inline expansion. 

Although inline expansion normally increases an executable file’s size, it can 
decrease executable size by avoiding the need to link in all or most of switched 
floating-point support 

With f 77, using x**y, mod, atan2, cabs, nint, and complex or doub- 
lecomplex arithmetic, can all provoke linking-in switched floating point even 
when -fsky, -f68881, or -ffpa is specified. Using fsky.il, 
f68881.il, or ffpa.il causes these calls to switched floating point to be 
replaced by inline code or calls to appropriate unswitched routines. 

With cc, use of almost any of the functions defined in <math . h> invokes 
switched floating point, but the appropriate . il file undoes the damage. With 
SP ICE 3 A . 7 , for instance, the version compiled with 
— — — — > 

-0 -ffpa /usr/lib/f fpa . il 

/ 
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was 40,000 bytes smaller than the version compiled just with 




Using the . il file also avoids certain overheads required for System V compa- 
tibility. 
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H 

System V Interface Compliance 



H.l. SVID History 



A joint goal of Sun and AT&T is to develop a version of UNIX that incorporates 
the best features of Berkeley BSD 4.2 and System V. Accordingly, Sun’s 
Software Release 3.2 mathematical library libm . a and related files have been 
modified so that C programs compiled with the -f sof t floating point option or 
compiled with -f switch and running without floating-point hardware better 
comply with the System V Interface Definition (SVID). Fortran and Pascal pro- 
grams, and C programs running with floating-point hardware, are usually not 
affected. The differences primarily involve exception handling. 

To understand the differences between exception handling according to SVID 
and the point of view represented by the IEEE Standard, it is necessary to review 
the circumstances under which both developed. Many of the ideas in SVID trace 
their origins to the early days of Unix, when it was first implemented on PDP- 
1 l’s and then ported to VAXes and IBM and Honeywell mainframe computers. 
These various environments have in common that rational floating-point opera- 
tions +, * and / are atomic machine instructions, while sqrt, conversion 

to integral value in floating-point format, and elementary transcendental func- 
tions are subroutines composed of many atomic machine instructions. 

Because these environments treat floating-point exceptions in varied ways, uni- 
formity could only be imposed by checking arguments and results in software 
before and after each atomic floating-point instruction. Since this would have too 
great an impact on performance, SVID does not specify the effect of floating- 
point exceptions such as division by zero or overflow. 

Operations implemented by subroutines are slow compared to single atomic 
floating-point instructions; extra error checking of arguments and results has little 
performance impact; so such checking is required by the SVID. When excep- 
tions are detected, default results are specified, errnoissetto EDOM for 
improper operands, or ERANGE for results that overflow or underflow, and the 
function matherr ( ) is called with a record containing details of the excep- 
tion. This costs little on the machines for which Unix was originally developed, 
but the value is correspondingly small since the far more common exceptions in 
the basic operations +, -, * and / are completely unspecified. 
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H.2. IEEE History 



H.3. SVID Future 
Directions 



H.4. Sun Implementation 



The IEEE Standard explicitly states that compatibility with previous implemen- 
tations was not a goal. Instead, an exception handling scheme was developed 
with efficiency and users’ requirements in mind. This scheme is uniform across 
the simple rational operations (+, * and /), and more complicated opera- 

tions such as remainder, square root, and conversion between formats. Although 
the Standard does not specify transcendental functions, the framers of the Stan- 
dard anticipated that the same exception handling scheme would be applied to 
elementary transcendental functions in conforming systems. 

Elements of IEEE exception handling include suitable default results and interr- 
uption of computation only when requested in advance. High performance con- 
forming processors such as the MC6888 1 compute the common elementary tran- 
scendental functions in single instructions that look, from a programmer’s 
viewpoint, just like the instructions for +, -, * and /. 

The current SVID identifies certain directions for future development. One of 
these is compatibility with the IEEE Standard. In particular a future version of 
the SVID will replace references to HUGE, intended to be a large finite number, 
with HUGE_VAL, which would be infinity on IEEE systems. HUGE_VAL 
would, for instance, be returned as the result of floating-point overflows. In this 
respect, Sun’s implementation has already arrived at SVID’s future direction, 
since the constant HUGE in /usr/include/math. h is defined to be a con- 
stant that compiles into IEEE infinity. 

In Release 3.2, the following C language libc and libm functions provide 
operand or result checking corresponding to SVID, when called from C programs 
compiled with -f soft, or compiled with -f switch and run without 
floating-point hardware: 

□ exp 

□ log and loglO 

□ pow 

□ sqrt 

□ hypot 
o cabs 

□ s inhand cosh 

□ sin, cos and tan 

□ asin, acos, atan and atan2 

When exceptional conditions are detected, the SVID function matherr ( ) is 
invoked. The default matherr () in libm returns a default result value and, 
for EDOM errors, prints an error message prior to continuing. Because SVID 
encompasses machines without infinities or NaNs, the default results specified 
are finite values and therefore sometimes misleading. Users may provide their 
own matherr ( ) function to obtain alternative processing. 

For efficiency, programs compiled with inline hardware floating-point or with 
-f switch and run with floating-point hardware do not do the extra checking 
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SIGNAL Notes 



libm . a Notes 



required to set EDOM or ERANGE or call matherr ( ) . Usually NaN is 
returned for the function value in situations where EDOM might otherwise be set, 
and Infinity is returned for the function value where ERANGE might be set to 
indicate overflow. Where ERANGE might be set to indicate underflow, these 
functions return subnormal numbers or zero. 

Note that if inline expansion files are used to expand libm functions, the SVID 
exception handling may be bypassed, even for -f sof t. 

The SVID defines SIGFPE as "floating-point exception". On Suns, SIGFPE is 
also generated by MC68010 and MC68020 integer division by zero, CHK and 
TRAPV instructions, and by FPA-related non-numerical exceptions. 

SVID specifies two floating-point exceptions, "PLOSS" (partial loss of 
significance) and "TLOSS" (total loss of significance). Unlike sqrt ( - 1 ) , 
these have no inherent mathematical meaning, and unlike exp(+-l0000), 
these do not reflect inherent limitations of a floating-point storage format. 

PLOSS and TLOSS reflect instead limitations of particular algorithms for f mod 
and for trigonometric functions that suffer abrupt declines in accuracy at definite 
boundaries. Like most IEEE implementations, the Sun algorithms used with 
-fsoft, -f68881,and -f f pa do not suffer such abrupt declines, and so do 
not signal PLOSS or TLOSS, nor do the Bessel functions which call the tri- 
gonometric functions. 

Instead, Sun’s sin, cos, and tan treat the essential singularity at infinity 
like other essential singularities by returning a NaN and setting EDOM for 
infinite arguments. The Bessel functions of the first kind, jO, jl, and jn, 
return zero for infinite arguments, while the Bessel functions of the second kind, 
yO, yl, and yn, return zero for positive infinite arguments, and return NaN and 
set EDOM for negative arguments. 

Likewise SVID specifies that fmod(x,y) should be zero if x/y would 
overflow, but Sun’s f mod, derived from the IEEE remainder function, does not 
compute x/y explicitly and hence always delivers a correcdy rounded result. 
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